
% Climate amenities paper model Q1_D2_FE0_S7_99.m
% This file runs the Original HetBin_150402.m first (PartOne) and then the bootstrap follows (PartTwo). 
% Base file completed on Thu, 2 April 2015
% Written by Hendrik Wolff and Ryan Kellogg

% When changing specifications, the following things need to change
 
% (1) save the Matlab .m file itself under the Specification of the regression,
% i.e. this one is Q1_D2_FE0_S7_99.m (then for Robustness, just add _R#
% depending on the number of the Robustness test, see Excel Table.)
% Within the code, change the file name for saving output using a find and replace. 
% Output file name shows up at the end of both part one and part two

% (2) Toggle selection changes in part one to whatever the new spec is.
% Part two for bootsrapping automatically uses the same toggles, though
% toggle_10 and toggle_11 are always set to 0


%------------------------------------------------------------

% Loading Toolbox for Hendrik
%cepath='M:\Private\Matlab\';
%path([cepath 'cetools;' cepath 'cedemos'],path);

% Loading toolbox for flux
%path('/home2/kelloggr/CompEcon/CEtools',path);
%path('/home2/kelloggr/CompEcon/CEdemos',path);

clear
format long
format compact

% diary Diary_Q1_D2_FE0_S7_99.txt

% Select main directory (comment out others)
% Ryan desktop
    HomeDir = 'C:\Work\ClimChgWelfare\Models\Models_2015_04';
% Hendrik desktop
    % HomeDir = 'R:\Project\Wolff\Research\ClimateChangeMichigan\dropboxFeb2012\Models_2012_02';
% Flux cluster
    % HomeDir = '/home2/kelloggr/ClimChgWelfare/Models_2015_04';

% Change to data directory
DirStr = [HomeDir '\Inputs_PUMA'];      % PC
% DirStr = [HomeDir '/Inputs_PUMA'];    % cluster
cd(DirStr)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD IN THE DATA AND GET SOME BASIC INFO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Puma identifiers and dependent variables (statefip PumaID msa 
% QOL QOL_orig QOL_25_1 QOL_GM QOL_Dahl QOL_fulladj QOL_doubleadj QOL_hs QOL_col Price Wage_orig Wage)
data = dlmread('QOL_ToMatlab.csv');  
K = size(data,2);       % # of columns
Stateraw = data(:,1); Pumaraw = data(:,2); msaraw = data(:,3);
DepVarraw = data(:,4:K);
clear data K
Nraw = size(Pumaraw,1);     % # of observations


% Climate normals: temperature bins (bin1 - bin222)
data = dlmread('ClimateBins_ToMatlab.csv');  
Binraw = data;
clear data
Nbin = size(Binraw,2);        % # of temperature bins


% Climate change 
data = dlmread('DeltaA2Ensemble2100_ToMatlab.csv');
DeltaA2Ensemble2100raw = data;
clear data
data = dlmread('DeltaA1F12100_ToMatlab.csv');
DeltaA1F12100raw = data;
clear data

% Geographic controls (slope SeaDistInv LakeDistInv)
data = dlmread('Controls_Geographic_ToMatlab.csv');
KCG0 = size(data,2);
XCGraw = data;
clear data

% Other weather (Precip RelHum Sun (all stations))
data = dlmread('Controls_OtherWx_RH_ToMatlab.csv');
KCOWA_RH = size(data,2);
XCOWA_RHraw = data;
clear data

% Other weather (Precip DewPt Sun (all stations))
data = dlmread('Controls_OtherWx_DP_ToMatlab.csv');
KCOWA_DP = size(data,2);
XCOWA_DPraw = data;
clear data

% Other weather 2 (Precip RelHumSummerTight Sun (all stations))
data = dlmread('Controls_OtherWx2_RH_ToMatlab.csv');
KCOWA2_RH = size(data,2);
XCOWA2_RHraw = data;
clear data

% Other weather 2 (Precip DewPtSummerTight Sun (all stations))
data = dlmread('Controls_OtherWx2_DP_ToMatlab.csv');
KCOWA2_DP = size(data,2);
XCOWA2_DPraw = data;
clear data

% Other weather 3 (PrecipSummer PrecipWinter RelHumSummer RelHumWinter SunSummer SunWinter (all stations))
data = dlmread('Controls_OtherWx3_RH_ToMatlab.csv');
KCOWA3_RH = size(data,2);
XCOWA3_RHraw = data;
clear data

% Other weather 3 (PrecipSummer PrecipWinter DewPtSummer DewPtWinter SunSummer SunWinter (all stations))
data = dlmread('Controls_OtherWx3_DP_ToMatlab.csv');
KCOWA3_DP = size(data,2);
XCOWA3_DPraw = data;
clear data

% Other weather (Precip RelHum Sun (4 stations))
data = dlmread('Controls_OtherWx4S_RH_ToMatlab.csv');
KCOW4S_RH = size(data,2);
XCOW4S_RHraw = data;
clear data

% Other weather (Precip DewPt Sun (4 stations))
data = dlmread('Controls_OtherWx4S_DP_ToMatlab.csv');
KCOW4S_DP = size(data,2);
XCOW4S_DPraw = data;
clear data

% Other weather 2 (Precip RelHumSummerTight Sun (4 stations))
data = dlmread('Controls_OtherWx4S2_RH_ToMatlab.csv');
KCOW4S2_RH = size(data,2);
XCOW4S2_RHraw = data;
clear data

% Other weather 2 (Precip DewPtSummerTight Sun (4 stations))
data = dlmread('Controls_OtherWx4S2_DP_ToMatlab.csv');
KCOW4S2_DP = size(data,2);
XCOW4S2_DPraw = data;
clear data

% Other weather 3 (PrecipSummer PrecipWinter RelHumSummer RelHumWinter SunSummer SunWinter (4 stations))
data = dlmread('Controls_OtherWx4S3_RH_ToMatlab.csv');
KCOW4S3_RH = size(data,2);
XCOW4S3_RHraw = data;
clear data

% Other weather 3 (PrecipSummer PrecipWinter DewPtSummer DewPtWinter SunSummer SunWinter (4 stations))
data = dlmread('Controls_OtherWx4S3_DP_ToMatlab.csv');
KCOW4S3_DP = size(data,2);
XCOW4S3_DPraw = data;
clear data

% Human / demographic controls (LogWeightedDensity Pct_HSEd Pct_BSEd Pct_GradEd Age PctHisp PctBlack)
data = dlmread('Controls_Demog_ToMatlab.csv');
KCH = size(data,2);
XCHraw = data;
clear data

% Unweighted population density (LogDensity)
data = dlmread('LogPopDens_ToMatlab.csv');
LogDenraw = data;
clear data

% Climate change projection for other weather (Precip RelHum Sun)
data = dlmread('DeltaA2Ensemble_OtherWx_RH_2100_ToMatlab.csv');
DeltaA2Ensemble2100_OW_RHraw = data;
clear data
data = dlmread('DeltaA1F1_OtherWx_RH_2100_ToMatlab.csv');
DeltaA1F12100_OW_RHraw = data;
clear data

% Climate change projection for other weather (Precip DewPt Sun)
data = dlmread('DeltaA2Ensemble_OtherWx_DP_2100_ToMatlab.csv');
DeltaA2Ensemble2100_OW_DPraw = data;
clear data
data = dlmread('DeltaA1F1_OtherWx_DP_2100_ToMatlab.csv');
DeltaA1F12100_OW_DPraw = data;
clear data

% Climate change projection for other weather 2 (Precip RelHumSummerTight Sun)
data = dlmread('DeltaA2Ensemble_OtherWx2_RH_2100_ToMatlab.csv');
DeltaA2Ensemble2100_OW2_RHraw = data;
clear data
data = dlmread('DeltaA1F1_OtherWx2_RH_2100_ToMatlab.csv');
DeltaA1F12100_OW2_RHraw = data;
clear data

% Climate change projection for other weather 2 (Precip DewPtSummerTight Sun)
data = dlmread('DeltaA2Ensemble_OtherWx2_DP_2100_ToMatlab.csv');
DeltaA2Ensemble2100_OW2_DPraw = data;
clear data
data = dlmread('DeltaA1F1_OtherWx2_DP_2100_ToMatlab.csv');
DeltaA1F12100_OW2_DPraw = data;
clear data

% Climate change projection for other weather 3 (PrecipSummer PrecipWinter RelHumSummer RelHumWinter SunSummer SunWinter)
data = dlmread('DeltaA2Ensemble_OtherWx3_RH_2100_ToMatlab.csv');
DeltaA2Ensemble2100_OW3_RHraw = data;
clear data
data = dlmread('DeltaA1F1_OtherWx3_RH_2100_ToMatlab.csv');
DeltaA1F12100_OW3_RHraw = data;
clear data

% Climate change projection for other weather 3 (PrecipSummer PrecipWinter DewPtSummer DewPtWinter SunSummer SunWinter)
data = dlmread('DeltaA2Ensemble_OtherWx3_DP_2100_ToMatlab.csv');
DeltaA2Ensemble2100_OW3_DPraw = data;
clear data
data = dlmread('DeltaA1F1_OtherWx3_DP_2100_ToMatlab.csv');
DeltaA1F12100_OW3_DPraw = data;
clear data

% load in income (for aggregating welfare impact)
data = dlmread('Income_ToMatlab.csv');
Incomeraw = data;
clear data

% load in population (for weighting) (Population Population_HS
% Population_Col)
data = dlmread('Population_ToMatlab.csv');
PopVecraw = data;
clear data

% load State FE
% below
data = dlmread('StateFE_ToMatlab.csv');
StateFEraw = data;
NFE = size(StateFEraw,2);
clear data

% load Division FE
% below
data = dlmread('DivisionFE_ToMatlab.csv');
DivFEraw = data;
NDFE = size(DivFEraw,2);
clear data

% load in bootstrapping group indicators (for clustered std errs)
data = dlmread('Bootstrap_ToMatlab.csv');
BSraw = data(:,2);      % 1st column for msa clustering, 2nd for state, 3rd for census division
clear data


% Go back to model directory
cd ..



%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOW SET UP ANALYSIS DATASET
% USE ALL LOCATIONS? INCLUDE CONTROLS? FIXED EFFECTS? 
% TOGGLES FOR THESE ARE IN THIS SECTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Toggle 1 is drop some msas: 0 = drops, 1 = no drops;
toggle_1 = 1;
    
% Toggle 2 is Dep Var: 1=QOL(commuting corrected) 2=QOL_orig, 3=QOL_25_1, 4=QOL_GM, 5=QOL_Dahl,
% 6=QOL_Notoadj, 7=QOL_doubleNotoadj, 8=QOL_hs, 9=QOL_col, 10=Price,
% 11=Wage_orig, 12=Wage (commuting corrected)
% for toggle_2 = 1:9
toggle_2 = 1;
    
% Toggle 3 is Geographic Controls: 0=None, 1=inv distance, 2=inv dist
% squared, 3=inv dist cubed, 4=inv dist squared + slope squared
% for toggle_3 = 1:4;
toggle_3 = 2;
    
% Toggle 4 is Other Weather Controls: 0=None, 1=Precip,
% 2=Precip+RelHum+Sunshine. 3=Precip+DewPt+Sunshine
% 4=Precip+RelHumSummer+Sunshine, 5=Precip+DewPtSummer+Sunshine, 6=Summer/Winter Precip+RelHum+Sunshine
% 7=Summer/Winter Precip+DewPt+Sunshine
% for toggle_4 = 0:7;
toggle_4 = 2;

% Toggle 4a is whether to use sunshine data interpolated from all stations
% or from the four nearest stations. 0 = all stations, 1 = four closest
% for toggle_4a = 0:1;
toggle_4a = 1;

% Toggle 5 is Human/demographic controls: 0=None, 1=Weighted density +
% education, 2=1+age+minorities, 3 = all vars squared
% for toggle_5 = 0:2:2;   
toggle_5 = 2;

% Toggle 6 is FE: 0=None, 1=Division, 2=State
toggle_6 = 0;

% Toggle 7 is Weights: 0=None, 1=Population
toggle_7 = 1;
    
% Toggle 8 is Kernel Weights: 0=heterog in temperature coefs only, 1=All
toggle_8 = 0;
        
% Toggle 9 is fraction of data to use for MWTP in cubic spline specification
toggle_9 = 0.99;
    
% Toggle 10 is whether to do set identification: 0=no, 1=yes (yes takes
% several hours to run)
toggle_10 = 0;

% Toggle 11 is whether to do the heterog pref standard errors and plots: 0 = no, 1 = yes
% Turn off if bootstrapping or looping over many toggles; this speeds things up
% No effect on point estimates.
toggle_11 = 1;

% Toggle 12: Number of spline basis functions. 3 means HDD / CDD. 
% 2 means HDD / CDD where valuations are equal
% 100 means 10 degree bins. 50 means 5 degree bins.
% 1 means 5th degree linear spline with kinks at 45 and 80 F
% 0 means 5th degree linear spline with kinks at 45 and 80 F; extremes are
% flat (so only 3df)
% -1 means 5th degree linear spline with kinks at 45 and 80 F; forced to be
% symmetric
toggle_12 = 7;

% Toggle 13: KERNEL BANDWIDTH
h = 2;


% Skip some combinations
%{
Skip = 0;
Skip = Skip + (toggle_1==0 & (toggle_2>1 | toggle_3~=1 | toggle_4~=3 | toggle_7==0 | toggle_8==0));
Skip = Skip + (toggle_2>2 & (toggle_1==0 | toggle_3~=1 | toggle_4~=3 | toggle_7==0 | toggle_8==0));
Skip = Skip + (toggle_3==0 & (toggle_5==1 | toggle_7==0 | toggle_8==0));
% Skip = Skip + (toggle_4<3 & (toggle_5==1 | toggle_7==0 | toggle_8==0));
Skip = Skip + (toggle_7==0 & (toggle_1==0 | toggle_2>2 | toggle_3~=1 | toggle_4~=3 | toggle_8==0));
Skip = Skip + (toggle_8==0 & (toggle_1==0 | toggle_2>2 | toggle_3~=1 | toggle_4~=3 | toggle_7==0));
if Skip>0
    continue;
end
%}

    
ID=cat(2, toggle_1, toggle_2, toggle_3, toggle_4, toggle_4a, toggle_5, toggle_6, toggle_7, toggle_8, toggle_9, toggle_10, toggle_11, toggle_12, h)
    
    % Toggle 4 indicates whether to use sunshine data interpolated from all
    % weather stations or just the 4 nearest stations
    if toggle_4a==0
        % Other weather (Precip RelHum Sun)
        KCOW_RH = KCOWA_RH; 
        XCOW_RHraw = XCOWA_RHraw;
        KCOW_DP = KCOWA_DP;
        XCOW_DPraw = XCOWA_DPraw; 
        % Other weather 2 (Precip RelHumSummerTight Sun)
        KCOW2_RH = KCOWA2_RH; 
        XCOW2_RHraw = XCOWA2_RHraw;
        KCOW2_DP = KCOWA2_DP;
        XCOW2_DPraw = XCOWA2_DPraw;
        % Other weather 3 (PrecipSummer PrecipWinter RelHumSummer RelHumWinter SunSummer SunWinter)
        KCOW3_RH = KCOWA3_RH; 
        XCOW3_RHraw = XCOWA3_RHraw;
        KCOW3_DP = KCOWA3_DP;
        XCOW3_DPraw = XCOWA3_DPraw;
    else
        % Other weather (Precip RelHum Sun)
        KCOW_RH = KCOW4S_RH; 
        XCOW_RHraw = XCOW4S_RHraw;
        KCOW_DP = KCOW4S_DP;
        XCOW_DPraw = XCOW4S_DPraw; 
        % Other weather 2 (Precip RelHumSummerTight Sun)
        KCOW2_RH = KCOW4S2_RH; 
        XCOW2_RHraw = XCOW4S2_RHraw;
        KCOW2_DP = KCOW4S2_DP;
        XCOW2_DPraw = XCOW4S2_DPraw;
        % Other weather 3 (PrecipSummer PrecipWinter RelHumSummer RelHumWinter SunSummer SunWinter)
        KCOW3_RH = KCOW4S3_RH; 
        XCOW3_RHraw = XCOW4S3_RHraw;
        KCOW3_DP = KCOW4S3_DP;
        XCOW3_DPraw = XCOW4S3_DPraw;
    end
    % delete stuff
    clear KCOW4S_RH XCOW4S_RHraw KCOW4S_DP XCOW4S_DPraw 
    clear KCOW4S2_RH XCOW4S2_RHraw KCOW4S2_DP XCOW4S2_DPraw 
    clear KCOW4S3_RH XCOW4S3_RHraw KCOW4S3_DP XCOW4S3_DPraw 
    clear KCOWA_RH XCOWA_RHraw KCOWA_DP XCOWA_DPraw 
    clear KCOWA2_RH XCOWA2_RHraw KCOWA2_DP XCOWA2_DPraw 
    clear KCOWA3_RH XCOWA3_RHraw KCOWA3_DP XCOWA3_DPraw    

    
    % Toggle 4 indicates whether to use dew point or relative humidity
    if toggle_4<=2 || toggle_4==4 || toggle_4==6
        KCOW = KCOW_RH; XCOWraw = XCOW_RHraw;
        KCOW2 = KCOW2_RH; XCOW2raw = XCOW2_RHraw;
        KCOW3 = KCOW3_RH; XCOW3raw = XCOW3_RHraw;
        DeltaA2Ensemble2100_OWraw = DeltaA2Ensemble2100_OW_RHraw;
        DeltaA1F12100_OWraw = DeltaA1F12100_OW_RHraw;
        DeltaA2Ensemble2100_OW2raw = DeltaA2Ensemble2100_OW2_RHraw;
        DeltaA1F12100_OW2raw = DeltaA1F12100_OW2_RHraw;
        DeltaA2Ensemble2100_OW3raw = DeltaA2Ensemble2100_OW3_RHraw;
        DeltaA1F12100_OW3raw = DeltaA1F12100_OW3_RHraw;
    else
        KCOW = KCOW_DP; XCOWraw = XCOW_DPraw;
        KCOW2 = KCOW2_DP; XCOW2raw = XCOW2_DPraw;  
        KCOW3 = KCOW3_DP; XCOW3raw = XCOW3_DPraw; 
        DeltaA2Ensemble2100_OWraw = DeltaA2Ensemble2100_OW_DPraw;
        DeltaA1F12100_OWraw = DeltaA1F12100_OW_DPraw;
        DeltaA2Ensemble2100_OW2raw = DeltaA2Ensemble2100_OW2_DPraw;
        DeltaA1F12100_OW2raw = DeltaA1F12100_OW2_DPraw;
        DeltaA2Ensemble2100_OW3raw = DeltaA2Ensemble2100_OW3_DPraw;
        DeltaA1F12100_OW3raw = DeltaA1F12100_OW3_DPraw;
    end
    clear KCOW_RH XCOW_RHraw KCOW2_RH XCOW2_RHraw KCOW3_RH XCOW3_RHraw
    clear DeltaA2Ensemble2100_OW_RHraw DeltaA1F12100_OW_RHraw
    clear DeltaA2Ensemble2100_OW2_RHraw DeltaA1F12100_OW2_RHraw
    clear DeltaA2Ensemble2100_OW3_RHraw DeltaA1F12100_OW3_RHraw
    clear KCOW_DP XCOW_DPraw KCOW2_DP XCOW2_DPraw KCOW3_DP XCOW3_DPraw
    clear DeltaA2Ensemble2100_OW_DPraw DeltaA1F12100_OW_DPraw 
    clear DeltaA2Ensemble2100_OW2_DPraw DeltaA1F12100_OW2_DPraw
    clear DeltaA2Ensemble2100_OW3_DPraw DeltaA1F12100_OW3_DPraw
    
    % Toggle 1 is dropping msas: 0 = drop some, 1 = no drops;
    if toggle_1==0
        Keep = find(msaraw~=160 & msaraw~=5602 & msaraw~=1122 & msaraw~=4472 & msaraw~=5520 & msaraw~=7362 & msaraw~=6162 & msaraw~=6320 & msaraw~=6480 & msaraw~=7120 & msaraw~=7480 & msaraw~=8000);
        %Keep = find(Stateraw~=6);
        USPersIncome2000=8576;
    elseif toggle_1==1
        Keep = [1:Nraw]';
        USPersIncome2000=8576;
    end
    
    % Now reduce dataset to just the observations corresponding to Keep
    N = length(Keep);
    State = Stateraw(Keep); Puma = Pumaraw(Keep); msa = msaraw(Keep); DepVar = DepVarraw(Keep,:);
    LogDen = LogDenraw(Keep);
    Bin = Binraw(Keep,:);
    DeltaA2Ensemble2100 = DeltaA2Ensemble2100raw(Keep,:);
    DeltaA1F12100 = DeltaA1F12100raw(Keep,:);
    XCG = XCGraw(Keep,:);
    XCOW = XCOWraw(Keep,:); XCOW2 = XCOW2raw(Keep,:); 
    XCOW3 = XCOW3raw(Keep,:); XCH = XCHraw(Keep,:);
    DeltaA2Ensemble2100_OW = DeltaA2Ensemble2100_OWraw(Keep,:);
    DeltaA2Ensemble2100_OW2 = DeltaA2Ensemble2100_OW2raw(Keep,:);
    DeltaA2Ensemble2100_OW3 = DeltaA2Ensemble2100_OW3raw(Keep,:);
    DeltaA1F12100_OW = DeltaA1F12100_OWraw(Keep,:);
    DeltaA1F12100_OW2 = DeltaA1F12100_OW2raw(Keep,:);
    DeltaA1F12100_OW3 = DeltaA1F12100_OW3raw(Keep,:);        
    IncomePC = Incomeraw(Keep);
    PopVec = PopVecraw(Keep,:);
    StateFE = StateFEraw(Keep,:); DivFE = DivFEraw(Keep,:);
    BS = BSraw(Keep);

    % get rid of FE for states that are no longer in the data
    blah = sum(StateFE);
    FEInd = find(blah~=0);
    StateFE = StateFE(:,FEInd);
    NFE = size(StateFE,2);
    clear Stateraw Pumaraw msaraw LogDenraw DepVarraw Binraw DeltaA2Ensemble2100raw DeltaA1F12100raw
    clear XCGraw XCOWraw XCOW2raw XCOW3raw XCHraw Incomeraw PopVecraw StateFEraw DivFEraw BSraw
    clear DeltaA2Ensemble2100_OWraw DeltaA2Ensemble2100_OW2raw DeltaA2Ensemble2100_OW3raw
    clear DeltaA1F12100_OWraw DeltaA1F12100_OW2raw DeltaA1F12100_OW3raw

    % Break out population data
    PopAll = PopVec(:,1);
    PopHS = PopVec(:,2);
    PopCol = PopVec(:,3);
    Pop = PopAll;           % Use total population for weighting. This will change if QOL_HS or QOL_Col are used
        
    % Puma-level total income
    Income = IncomePC .* Pop;   % Use total income for weighting. This will change if QOL_HS or QOL_Col are used
    

    % DROP ALL TEMPERATURE BINS THAT ARE NEVER IN THE DATA (EXTREME COLD AND
    % EXTREME HOT DAYS)
    temp = max(Bin) + max(DeltaA2Ensemble2100) + abs(min(DeltaA2Ensemble2100));
    temp = temp + max(DeltaA1F12100) + abs(min(DeltaA1F12100));     % will not be zero if any activity
    Ind = find(temp>0);        % relevant range of data
    MinInd = min(Ind'); MaxInd = max(Ind');
    TminC = -50 + 0.5*(MinInd-2) + 0.25;  % min temp in relevant range
    TmaxC = -50 + 0.5*(MaxInd-2) + 0.25;  % max temp in relevant range
    Tmin = TminC * 1.8 + 32;
    Tmax = TmaxC * 1.8 + 32;
    TempsC = [TminC:0.5:TmaxC]';     % vector of temperatures
    Temps = TempsC * 1.8 + 32;      % in degrees F
    Inds = [MinInd:1:MaxInd];       % full covered range (no holes)
    Bin = Bin(:,Inds); DeltaA2Ensemble2100 = DeltaA2Ensemble2100(:,Inds); 
    DeltaA1F12100 = DeltaA1F12100(:,Inds); 
    Nbin = size(Bin,2);
    clear temp Ind MinInd MaxInd Inds

    % CALCULATE FUTURE TEMPERATURES
    BinA2Ensemble2100 = Bin + DeltaA2Ensemble2100;
    BinA1F12100 = Bin + DeltaA1F12100;


    % GET AVERAGE US TEMPERATURES (POPULATION WEIGHTED)
    TempAvgPres = (Bin * Temps)' * Pop / sum(Pop) / 365;
    TempAvgA2Ensemble = (BinA2Ensemble2100 * Temps)' * Pop / sum(Pop) / 365;
    TempAvgA1F1 = (BinA1F12100 * Temps)' * Pop / sum(Pop) / 365;

    % Identify some specific temperature bins
    Bin0 = floor((0-Tmin+0.449)/0.9) + 1;
    Bin10 = floor((10-Tmin+0.449)/0.9) + 1; % adding 0.449 means that a temp on borderline between bins goes to lower bin
    Bin20 = floor((20-Tmin+0.449)/0.9) + 1; % 50F is on borderline, this forces it to lower bin. Makes 50F same number of bins away from 65 as is 80
    Bin25 = floor((25-Tmin+0.449)/0.9) + 1;
    Bin30 = floor((30-Tmin+0.449)/0.9) + 1;
    Bin35 = floor((35-Tmin+0.449)/0.9) + 1;
    Bin40 = floor((40-Tmin+0.449)/0.9) + 1; 
    Bin45 = floor((45-Tmin+0.449)/0.9) + 1; 
    Bin50 = floor((50-Tmin+0.449)/0.9) + 1; 
    Bin55 = floor((55-Tmin+0.449)/0.9) + 1;     
    Bin60 = floor((60-Tmin+0.449)/0.9) + 1;     
    Bin65 = floor((65-Tmin+0.449)/0.9) + 1; 
    Bin70 = floor((70-Tmin+0.449)/0.9) + 1; 
    Bin75 = floor((75-Tmin+0.449)/0.9) + 1; 
    Bin80 = floor((80-Tmin+0.449)/0.9) + 1; 
    Bin85 = floor((85-Tmin+0.449)/0.9) + 1; 
    Bin90 = floor((90-Tmin+0.449)/0.9) + 1; 
    Bin95 = floor((95-Tmin+0.449)/0.9) + 1; 
    Bin100 = floor((100-Tmin+0.449)/0.9) + 1; 
    Bin110 = floor((110-Tmin+0.449)/0.9) + 1; 
    
    % FORM SPLINES / DUMMIES / HDD & CDD ON THE TEMPERATURE BINS
    % NUMBER OF SPLINE BASIS FUNCTIONS IS S
    % THIS SECTION CONCLUDES WITH S COVARIATES IN THE X MATRIX
    % USES THE MIRANDA-FACKLER TOOLBOX
    % NOTE THAT SPLINES DO NOT REQUIRE A CONSTANT IN THE REGRESSION!

    %fspace1 = fundef({'spli',[-60.5 0 40 60 70 80 140]});
    %xnode1 = funnode(fspace1)
    %fspace1 = fundef({'spli',[-60.5 140.5]})

    if toggle_12==50           % 5 degree temperature bins
        S = 20;     % 20 5-degree bins. <5 and >=95 grouped
        Basis = zeros(Nbin,S);      % basis function for binned preferences
        for s = 1:S         % Loop to fill in the basis matrix
            if s==1         % T<5
                ind = find(Temps<5);
                Basis(ind,s) = 1;
            elseif s==20        %T>=95
                ind = find(Temps>=95);
                Basis(ind,s) = 1;
            else
                bot = s * 5 - 5;
                top = bot + 5;
                ind = find(Temps>=bot & Temps<top);
                Basis(ind,s) = 1;
            end
        end
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates     
        Pctile = 1;                 % flag that Basis covers all temperatures
    elseif toggle_12==100           % 10 degree temperature bins
        S = 7;     % 10 10-degree bins. <35 and >=85 grouped
        Basis = zeros(Nbin,S);      % basis function for binned preferences
        for s = 1:S         % Loop to fill in the basis matrix
            if s==1         % T<35
                ind = find(Temps<35);
                Basis(ind,s) = 1;
            elseif s==7        %T>=85
                ind = find(Temps>=85);
                Basis(ind,s) = 1;
            else
                bot = s * 10 + 15;
                top = bot + 10;
                ind = find(Temps>=bot & Temps<top);
                Basis(ind,s) = 1;
            end
        end
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates  
        Pctile = 1;                 % flag that Basis covers all temperatures
    elseif toggle_12>=4 & toggle_12<50        % use Miranda-Fackler splines
        S = toggle_12;
        % Define splines only over the range of present-day temperatures
        temp = max(Bin);
        Ind = find(temp>0); 
        MinInd = min(Ind'); MaxInd = max(Ind');
        Inds = [MinInd:1:MaxInd];       % full covered range (no holes)
        BinInd = Bin(:,Inds);           % bins spanning present-day data
        TempsInd = Temps(Inds);
        TminInd = min(TempsInd); TmaxInd = max(TempsInd);
        NbinInd = length(Inds);         % number of 0.5C bins in selected range  
        % Find population-weighted temperature distribution
        BinPop = Bin .* repmat(Pop,1,Nbin);
        DaycountPopWt = sum(BinPop) ./ repmat(sum(Pop),1,Nbin);
        CumDays = zeros(1,Nbin);    % initialize: accumulates days over bins
        CumDays(1) = DaycountPopWt(1);
        for i = 2:Nbin;
            CumDays(i) = CumDays(i-1) + DaycountPopWt(i);    % accumulate
        end
        TotDays = max(CumDays);
        CumDays = CumDays / TotDays;                % cdf
        % Generate even spacing on percentiles of the temperature distribution
        Step = 1 / (S - 3);
        BreaksPct = 0:Step:1;       % breakpoints in percentiles
        BreaksTemp = zeros(1,S-2);  % initialize breakpoints in temperatures
        BreaksTemp(1) = TminInd; BreaksTemp(S-2) = TmaxInd; % first and last breaks
        if S>4
            for b = 2:S-3
                Ind = min([find(CumDays>=BreaksPct(b))]');
                BreaksTemp(b) = Temps(Ind);
            end
        end
        %fspace = fundefn('spli',S,TminInd,TmaxInd);       % set up basis functions
        fspace = fundef({'spli',BreaksTemp});
        Basis = funbas(fspace,TempsInd);       % basis functions: NbinInd * S matrix              
        X = [BinInd * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates
        % Define 5th order spline for the local linear pref weighting
        % matrix. For use if S>5 & S<50
        SW = 5;
        StepW = 1 / (SW - 3);
        BreaksPctW = 0:StepW:1;
        BreaksTempW = zeros(1,SW-2);
        BreaksTempW(1) = TminInd; BreaksTempW(SW-2) = TmaxInd; % first and last breaks
        if SW>4
            for b = 2:SW-3
                Ind = min([find(CumDays>=BreaksPctW(b))]');
                BreaksTempW(b) = Temps(Ind);
            end
        end
        fspaceW = fundef({'spli',BreaksTempW});
        BasisW = funbas(fspaceW,TempsInd);       % basis functions: NbinInd * S matrix              
        XW = [BinInd * BasisW];          % covariates: N * S matrix
    elseif toggle_12==3
        S = 3;       % HDD / CDD regressions (with constant)
        Break1 = Bin65;  % index of bin containing max MWTP
        Basis = zeros(Nbin,3);  % initialize
        Basis(:,1) = 1;         % constant
        for i = 1:Break1
            Basis(i,2) = Break1 - i;     % line on HDD side
        end
        for i = Break1:Nbin
            Basis(i,3) = i - Break1;     % line on CDD side
        end
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates
        Pctile = 1;                 % flag that Basis covers all temperatures
    elseif toggle_12==2
        S = 2;       % HDD / CDD regressions (with constant, B_HDD = B_CDD)
        Break1 = Bin65;  % index of bin containing max MWTP
        Basis = zeros(Nbin,2);  % initialize
        Basis(:,1) = 1;         % constant
        for i = 1:Break1
            Basis(i,2) = Break1 - i;     % line on HDD side
        end
        for i = Break1:Nbin
            Basis(i,2) = i - Break1;     % line on CDD side
        end
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates
        Pctile = 1;                 % flag that Basis covers all temperatures
    elseif toggle_12==1
        S = 5;      % 5th degree linear spline
        Break1 = Bin45;
        Break2 = Bin65;
        Break3 = Bin80;
        Basis = zeros(Nbin,5);  % initialize
        Basis(:,1) = 1;         % constant
        % Build bases
        for i = 1:Break1
            Basis(i,2) = i; 
        end
        for i = Break1+1:Break2
            Basis(i,2) = Break1;
            Basis(i,3) = i-Break1;
        end
        for i = Break2+1:Break3
            Basis(i,2) = Break1;
            Basis(i,3) = Break2-Break1;
            Basis(i,4) = i-Break2;
        end
        for i = Break3+1:Nbin
            Basis(i,2) = Break1;
            Basis(i,3) = Break2-Break1;
            Basis(i,4) = Break3-Break2;   
            Basis(i,5) = i-Break3;
        end
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates
        Pctile = 1;                 % flag that Basis covers all temperatures  
    elseif toggle_12==0
        S = 3;      % 3rd degree linear spline--extremes forced flat
        Break1 = Bin45;
        Break2 = Bin65;
        Break3 = Bin80;
        Basis = zeros(Nbin,3);  % initialize
        Basis(:,1) = 1;         % constant
        % Build bases
        for i = Break1+1:Break2
            Basis(i,2) = i-Break1;
        end
        for i = Break2+1:Break3
            Basis(i,2) = Break2-Break1;
            Basis(i,3) = i-Break2;
        end
        for i = Break3+1:Nbin
            Basis(i,2) = Break2-Break1;
            Basis(i,3) = Break3-Break2;
        end
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates
        Pctile = 1;                 % flag that Basis covers all temperatures     
    elseif toggle_12==-1
        S = 3;      % 5th degree linear spline--forced symmetric
        Break1 = Bin45;
        Break2 = Bin65;
        Break3 = Bin80;
        Basis = zeros(Nbin,3);  % initialize
        Basis(:,1) = 1;         % constant   
        % Build bases
        for i = 1:Break1
            Basis(i,2) = i; 
        end
        for i = Break1+1:Break2
            Basis(i,2) = Break1;
            Basis(i,3) = i-Break1;
        end
        for i = Break2+1:Break3
            Basis(i,2) = Break1;
            Basis(i,3) = Break2-Break1-(i-Break2);
        end
        for i = Break3+1:Nbin
            Basis(i,2) = Break1-(i-Break3);
            Basis(i,3) = Break2-Break1-(Break3-Break2);
        end       
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates
        Pctile = 1;                 % flag that Basis covers all temperatures  
    elseif toggle_12==-2
        S = 5;      % 5th degree linear spline--extremes forced flat.
        Break1 = Bin25;
        Break1 = floor((4.55-Tmin+0.449)/0.9) + 1;
        Break2 = Bin55;
        Break3 = Bin65;   
        Break4 = Bin80; 
        Break5 = Bin90; 
        Break5 = floor((87.35-Tmin+0.449)/0.9) + 1;
        Basis = zeros(Nbin,5);  % initialize
        Basis(:,1) = 1;         % constant
        % Build bases
        for i = Break1+1:Break2
            Basis(i,2) = i-Break1;
        end        
        for i = Break2+1:Break3
            Basis(i,2) = Break2-Break1;
            Basis(i,3) = i-Break2;
        end
        for i = Break3+1:Break4
            Basis(i,2) = Break2-Break1;
            Basis(i,3) = Break3-Break2;
            Basis(i,4) = i-Break3;
        end
        for i = Break4+1:Break5
            Basis(i,2) = Break2-Break1;
            Basis(i,3) = Break3-Break2;     
            Basis(i,4) = Break4-Break3;
            Basis(i,5) = i-Break4;
        end
        for i = Break5+1:Nbin
            Basis(i,2) = Break2-Break1;
            Basis(i,3) = Break3-Break2;     
            Basis(i,4) = Break4-Break3;
            Basis(i,5) = Break5-Break4;
        end        
        X = [Bin * Basis];          % covariates: N * S matrix
        K = S;                      % # of covariates
        Pctile = 1;                 % flag that Basis covers all temperatures
    end
    WCols = ones(1,K);      % initialize vector of which columns to use in weighted regressions

    % Clarify number of bins used in estimation and location of 65F
    if toggle_12<=3 || toggle_12>=50        % not using cubic splines; basis covers all temperatures
        NbinInd = Nbin;     % basis covers all temperatures
        Bin65Ind = Bin65;
        Bin40Ind = Bin40;
        Bin80Ind = Bin80;
    else                % basis covers short range
        Bin65Ind = floor((65-TminInd+0.449)/0.9) + 1;    % index of the bin containing 65F
        Bin40Ind = floor((40-TminInd+0.449)/0.9) + 1;    % index of the bin containing 40F
        Bin80Ind = floor((80-TminInd+0.449)/0.9) + 1;    % index of the bin containing 80F
    end
    
   
    
    % Toggle 2 is Dep Var: 1=QOL(commuting corrected) 2=QOL_orig, 3=QOL_25_1, 4=QOL_GM, 5=QOL_Dahl,
    % 6=QOL_Notoadj, 7=QOL_doubleNotoadj, 8=QOL_hs, 9=QOL_col, 10=Price,
    % 11=Wage_orig, 12=Wage (commuting corrected)
    if toggle_2==1
        Q = DepVar(:,1);     % QOL, with wages at place of work and commuting adjusted
    elseif toggle_2==2
        Q = DepVar(:,2);     % QOL_orig
    elseif toggle_2==3
        Q = DepVar(:,3);     % Old QOL: 0.25p - 2
    elseif toggle_2==4
        Q = DepVar(:,4);     % Glaeser-Mare correction: 1/3 wage diff is selection
    elseif toggle_2==5
        Q = DepVar(:,5);     % QOL adjusted for migration/selection per Dahl(2002)
    elseif toggle_2==6
        Q = DepVar(:,6);     % QOL adjusted for migration per Noto
    elseif toggle_2==7
        Q = DepVar(:,7);     % QOL double adjusted for migration per Noto
    elseif toggle_2==8
        Q = DepVar(:,8);     % QOL for high school education people
        Pop = PopNonCol;
        Income = Income .* PopNonCol ./ PopAll;     % scale income
    elseif toggle_2==9
        Q = DepVar(:,9);     % QOL for college education people
        Pop = PopCol;
        Income = Income .* PopCol ./ PopAll;     % scale income
    elseif toggle_2==10
        Q = DepVar(:,10);     % house prices
    elseif toggle_2==11
        Q = -DepVar(:,11);    % original wages
    elseif toggle_2==12
        Q = -DepVar(:,12);    % wages adjusted for commuting
    end

    % Toggle 3 is Geographic Controls: 0=None, 1=inv distance, 2=inv dist
    % squared, 3=inv dist cubed, 4=inv dist squared + slope squared
    if toggle_3==0
        KCG = 0;        % no controls
    elseif toggle_3==1  % Keep Great Lakes and higher order sea distances out of kernel weights
        WCols = [WCols 1 1 0];
        X = [X XCG(:,1:3)]; K = size(X,2); KCG = size(XCG(:,1:3),2);
    elseif toggle_3==2
        WCols = [WCols 1 1 0 0 0];
        X = [X XCG(:,1:5)]; K = size(X,2); KCG = size(XCG(:,1:5),2);
    elseif toggle_3==3
        WCols = [WCols 1 1 0 0 0 0 0];
        X = [X XCG(:,1:7)]; K = size(X,2); KCG = size(XCG(:,1:7),2);
    elseif toggle_3==4
        WCols = [WCols 1 1 0 0 0 0];
        X = [X XCG(:,[1:5 8])]; K = size(X,2); KCG = size(XCG(:,[1:5 8]),2);
    end
    
    % Toggle 4 is Other Weather Controls: 0=None, 1=Precip,
    % 2=Precip+RelHum+Sunshine. 3=Precip+DewPt+Sunshine
    % 4=Precip+RelHumSummer+Sunshine, 5=Precip+DewPtSummer+Sunshine, 6=Summer/Winter Precip+RelHum+Sunshine
    % 7=Summer/Winter Precip+DewPt+Sunshine
    if toggle_4==0
    elseif toggle_4==1
        X = [X XCOW(:,1)]; K = size(X,2);
        WCols = [WCols 1];
    elseif toggle_4==2 || toggle_4==3
        X = [X XCOW]; K = size(X,2);
        WCols = [WCols 1 1 1];
    elseif toggle_4==4 || toggle_4==5
        X = [X XCOW2]; K = size(X,2);
        WCols = [WCols 1 1 1];     
    elseif toggle_4>=6
        X = [X XCOW3]; K = size(X,2);
        WCols = [WCols 1 1 1 1 1 1];
    end
    
    % Toggle 5 is Human/demographic controls: 0=None, 1=Some, 2=All
    if toggle_5==0
    elseif toggle_5==1
        X = [X XCH(:,1:4)]; K = size(X,2);
        WCols = [WCols 1 1 1 1];
    elseif toggle_5==2
        X = [X XCH]; K = size(X,2);  
        WCols = [WCols ones(1,KCH)];
    elseif toggle_5==3
        X = [X XCH XCH.^2]; K = size(X,2);  
        WCols = [WCols ones(1,KCH*2)];
    end
    
    % Toggle 6 is FE: 0=None, 1=Division, 2=State
    if toggle_6==0
        FE = 0;
        DFE = 0;
    elseif toggle_6==1
        FE = 0;
        DFE = 1;
    elseif toggle_6==2
        FE = 1;
        DFE = 0;
    end

    % Toggle 7 is Weights: 0=None, 1=Population
    PopStored = Pop;
    if toggle_7==0
        Pop = ones(N,1);
    elseif toggle_7==1
        Pop = Pop / 100000;     % normalize
    end
    
    % Toggle 8 is Kernel Weights: 0=heterog prefs on temps only, 1=All
    blah = toggle_3 + toggle_4 + toggle_5;
    if toggle_8==0 && blah>0
        WCols(S+1:K) = 0;
    else
    end
       
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % HOMOGENOUS PREFERENCE REGRESSION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if FE==0 && DFE==0
        XAll = X;               % no FE.  
    elseif FE==0 && DFE==1
        XAll = [X DivFE];       % add division FE
    else
        XAll = [X StateFE];     % add state FE
    end
    % Population weighing
    XXX = diag(Pop.^0.5) * XAll;
    QQQ = diag(Pop.^0.5) * Q;
    
    KKK = size(XXX,2);        % # of vars in regression
    Qinv = inv(XXX' * XXX);
    BetaHomog = Qinv * XXX' * QQQ;
        
    % Variance of BetaHomog
    Resids = QQQ - XXX * BetaHomog;
    XXXe = XXX .* repmat(Resids,1,KKK);
    VarBetaHomogNoClust = N / (N - KKK) * (Qinv * XXXe' * XXXe * Qinv); % robust variance
    % Now do clustering: Qinv * X'ee'X * Qinv. ee' has non-zero off-diagonals
    % for common msas. Start by calculating this ee' matrix (called EE)
    EE = zeros(N,N); 
    % EE = sparse(EE);
    for r = 1:N         % loop over rows of EE
        target = BS(r);     % msa of row r
        Ind = find(BS==target);     % all columns with msa equal to target
        EE(r,Ind) = 1;
        EE(r,:) = Resids(r) * Resids' .* EE(r,:);   % only capture residuals common to msa
    end
    VarBetaHomog = N / (N - KKK) * (Qinv * XXX' * EE * XXX * Qinv); % clustered variance
    clear EE        % it's big, get rid of it
    % VarBetaHomog = VarBetaHomogNoClust;  UNCOMMENT TO USE NON-CLUSTERED VARIANCE
    SDBetaHomog = (diag(VarBetaHomog)).^0.5;
    MeanQQQ = mean(QQQ); MeanResids = mean(Resids);
    R2 = 1 - (Resids-MeanResids)'*(Resids-MeanResids) / ((QQQ-MeanQQQ)'*(QQQ-MeanQQQ));
    SSE = (Resids-MeanResids)'*(Resids-MeanResids);
    
    
    % Get value of an extra day at 65 degrees. Will need to subtract this off
    % of all estimates
    Val65h = Basis(Bin65Ind,:) * BetaHomog(1:S);

    % Marginal impact
    BinValueh = Basis * BetaHomog(1:S) - repmat(Val65h,NbinInd,1);

    % Std error of marginal impact
    BasisDiff = Basis - repmat(Basis(Bin65Ind,:),NbinInd,1);
    VarBinValueh = zeros(NbinInd,1);       % initialize
    for b = 1:NbinInd
        VarBinValueh(b) = BasisDiff(b,:) * VarBetaHomog(1:S,1:S) * BasisDiff(b,:)';
    end
    SDBinValueh = VarBinValueh.^0.5;
    TopInth = BinValueh + 1.96 * SDBinValueh;   % top of 95% c.i.
    BotInth = BinValueh - 1.96 * SDBinValueh;   % bottom of 95% c.i.
    
    
    % Population-weighted temperature distributions (for plotting)
    BinPop = Bin .* repmat(PopStored,1,Nbin);
    Daycounth = sum(BinPop) ./ repmat(sum(PopStored),1,Nbin);
    BinA2Ensemble2100Pop = BinA2Ensemble2100 .* repmat(PopStored,1,Nbin);
    DaycountA2Ensemble2100h = sum(BinA2Ensemble2100Pop) ./ repmat(sum(PopStored),1,Nbin);
    BinA1F12100Pop = BinA1F12100 .* repmat(PopStored,1,Nbin);
    DaycountA1F12100h = sum(BinA1F12100Pop) ./ repmat(sum(PopStored),1,Nbin);    
        
    % If using splines, extrapolate MWTP at extremes
    if toggle_12>3 && toggle_12<50
        % Initialize valuations; BinValueh within present-day temp range, zero outside
        BinValuehInd = BinValueh;     % defined over present-day temp range
        BinValueh = zeros(Nbin,1);    % defined over present & future temps (wider)
        BinValueh(Inds') = BinValuehInd;
        % do conf interval & std errors too
        TopInthInd = TopInth; BotInthInd = BotInth;
        TopInth = zeros(Nbin,1); BotInth = zeros(Nbin,1);
        TopInth(Inds') = TopInthInd; BotInth(Inds') = BotInthInd;
        SDBinValuehInd = SDBinValueh;
        SDBinValueh = zeros(Nbin,1);
        SDBinValueh(Inds') = SDBinValuehInd;
        % Define range of temperatures over which splines are still part of MWTP
        Pctile = 1 - (1 - toggle_9) / 2;    % percent of distribution (income-weighted) to exclude on each side
        CenterSh = find(CumDays>(1-Pctile) & CumDays<Pctile);    % select center of distribution
        LowSh = find(CumDays<=(1-Pctile)); HighSh = find(CumDays>=Pctile);
        MinSh = min(CenterSh'); MaxSh = max(CenterSh'); 
        % Now flatten valuations outside of the active spline range
        ColdVal = BinValueh(MinSh); HotVal = BinValueh(MaxSh);    % valuations at bottom and top of temp dist
        BinValueh(LowSh) = ColdVal; BinValueh(HighSh) = HotVal;
        % do conf interval too
        TopInth(LowSh) = TopInth(MinSh); TopInth(HighSh) = TopInth(MaxSh);
        BotInth(LowSh) = BotInth(MinSh); BotInth(HighSh) = BotInth(MaxSh);
    end

    
    % Plot marginal impact from homogenous regression
    tempminplot = Temps(Bin0); tempmaxplot = Temps(Bin110);
    temp = [tempminplot:0.9:tempmaxplot];  % vector of temperatures to plot
    clf
    [AX,H1,H2] = plotyy(temp,BinValueh(Bin0:Bin110),temp,Daycounth(Bin0:Bin110));
    set(AX(1),'YLim',[-0.002 0.0005]); set(AX(2),'YLim',[0 20]);
    y1inc = (0.002 + 0.0005) / 5; y2inc = (20 - 0) / 5;
    set(AX(1),'YTick',[-0.002:y1inc:0.0005]); set(AX(2),'YTick',[0:y2inc:20]);
    grid on
    set(H1,'Color','b'); set(H1,'LineWidth',3);      % preference estimate
    set(H2,'LineStyle','--');  set(H2,'LineWidth',2); set(H2,'Color','k');  % current temps
    axes(AX(1)); ylabel('Marginal Value of an Extra Day at Each Temperature');
    axes(AX(2)); ylabel('Pop-weighted average number of Days at Each Temperature');
    % add 2050 temps, 2100 temps, top of c.i., and bottom of c.i.
    H3 = line(temp,DaycountA2Ensemble2100h(Bin0:Bin110),'LineStyle','-.','LineWidth',2,'Color','r','Parent',AX(2));
    % H3 = line(temp,DaycountA1F12100h(Bin0:Bin110),'LineStyle','-.','LineWidth',2,'Color','g','Parent',AX(2));
    H4 = line(temp,TopInth(Bin0:Bin110),'LineStyle','--','LineWidth',2,'Color','b','Parent',AX(1));
    H5 = line(temp,BotInth(Bin0:Bin110),'LineStyle','--','LineWidth',2,'Color','b','Parent',AX(1));
    set(AX(1),'XLim',[0 110]); set(AX(1),'XTick',[0:10:110]);
    set(AX(2),'XLim',[0 110]); set(AX(2),'XTick',[0:10:110]);
   
    
    
    % TEST OF H0: BETA_HDD = BETA_CDD
    if toggle_12==3
        TestStat = abs(BetaHomog(2)-BetaHomog(3)) / ([1 -1] * VarBetaHomog(2:3,2:3) * [1; -1])^0.5;
        TestPVal = 2 * (1-normcdf(TestStat,0,1));
    elseif toggle_12==0
        TestStat = abs(BetaHomog(2)+BetaHomog(3)) / ([1 1] * VarBetaHomog(2:3,2:3) * [1; 1])^0.5;
        TestPVal = 2 * (1-normcdf(TestStat,0,1));
    elseif toggle_12==1
        TestStat = abs(BetaHomog(3)+BetaHomog(4)) / ([1 1] * VarBetaHomog(3:4,3:4) * [1; 1])^0.5;
        TestPVal = 2 * (1-normcdf(TestStat,0,1));
    end
    
    % TEST OF H0: MODEL IS HDD/CDD, NOT 5 PART LINEAR SPLINE
    if toggle_12==1
        R = [0 1 -1 0 0; 0 0 0 1 -1]; J = 2;
        TestStat3 = (R*BetaHomog(1:5))' * inv(R * VarBetaHomog(1:5,1:5) * R') * R*BetaHomog(1:5) / J;
        DF = N - K;
        TestPVal3 = 1 - fcdf(TestStat3,J,DF);
        R = [0 1 -1 0 0]; J = 1;
        TestStat3Cold = (R*BetaHomog(1:5))' * inv(R * VarBetaHomog(1:5,1:5) * R') * R*BetaHomog(1:5) / J;
        DF = N - K;
        TestPVal3Cold = 1 - fcdf(TestStat3Cold,J,DF);
        R = [0 0 0 1 -1]; J = 1;
        TestStat3Heat = (R*BetaHomog(1:5))' * inv(R * VarBetaHomog(1:5,1:5) * R') * R*BetaHomog(1:5) / J;
        DF = N - K;
        TestPVal3Heat = 1 - fcdf(TestStat3Heat,J,DF);
    end
    % TEST OF H0: MODEL IS FLAT AT EXTREMES, NOT 5 PART LINEAR SPLINE
    if toggle_12==1
        R = [0 1 0 0 0; 0 0 0 0 1]; J = 2;
        TestStat0 = (R*BetaHomog(1:5))' * inv(R * VarBetaHomog(1:5,1:5) * R') * R*BetaHomog(1:5) / J;
        DF = N - K;
        TestPVal0 = 1 - fcdf(TestStat0,J,DF);
        R = [0 1 0 0 0]; J = 1;
        TestStat0Cold = (R*BetaHomog(1:5))' * inv(R * VarBetaHomog(1:5,1:5) * R') * R*BetaHomog(1:5) / J;
        DF = N - K;
        TestPVal0Cold = 1 - fcdf(TestStat0Cold,J,DF);
        R = [0 0 0 0 1]; J = 1;
        TestStat0Heat = (R*BetaHomog(1:5))' * inv(R * VarBetaHomog(1:5,1:5) * R') * R*BetaHomog(1:5) / J;
        DF = N - K;
        TestPVal0Heat = 1 - fcdf(TestStat0Heat,J,DF);
    end    
    
    
    
    %%%%%%%%
    % HOMOGENOUS PREF WELFARE CALCULATIONS
    % Start with welfare impact of temperature change
    DQA2Ensemble2100h_T = DeltaA2Ensemble2100 *  BinValueh;   % vector of delta welfare as pct of income
    DQA1F12100h_T = DeltaA1F12100 *  BinValueh;   % vector of delta welfare as pct of income

    AggWelfareA2Ensemble2100h_T = DQA2Ensemble2100h_T' * Income;          % in billions
    AggWelfareA1F12100h_T = DQA1F12100h_T' * Income;          % in billions
    % as a percent of total income
    WelfarePctA2Ensemble2100h_T = AggWelfareA2Ensemble2100h_T / sum(Income);
    WelfarePctA1F12100h_T = AggWelfareA1F12100h_T / sum(Income);
    % Get welfare impact in $ using total US personal income
    AggWelfareA2Ensemble2100h_T = WelfarePctA2Ensemble2100h_T * USPersIncome2000;
    AggWelfareA1F12100h_T = WelfarePctA1F12100h_T * USPersIncome2000;
    
    % Break into winter gain vs. summer loss
    DQA2Ensemble2100h_TW = DeltaA2Ensemble2100(:,1:Bin65) *  BinValueh(1:Bin65);   % vector of delta welfare as pct of income
    DQA1F12100h_TW = DeltaA1F12100(:,1:Bin65) *  BinValueh(1:Bin65);   % vector of delta welfare as pct of income    
    DQA2Ensemble2100h_TS = DeltaA2Ensemble2100(:,Bin65+1:Nbin) *  BinValueh(Bin65+1:Nbin);   % vector of delta welfare as pct of income
    DQA1F12100h_TS = DeltaA1F12100(:,Bin65+1:Nbin) *  BinValueh(Bin65+1:Nbin);   % vector of delta welfare as pct of income    
    WelfarePctA2Ensemble2100h_TW = DQA2Ensemble2100h_TW' * Income / sum(Income);
    WelfarePctA2Ensemble2100h_TS = DQA2Ensemble2100h_TS' * Income / sum(Income);
    WelfarePctA1F12100h_TW = DQA1F12100h_TW' * Income / sum(Income);
    WelfarePctA1F12100h_TS = DQA1F12100h_TS' * Income / sum(Income);
    
    % Welfare impact of precip change
    if toggle_4<1       % precipitation not in model
        DQA2Ensemble2100h_P = zeros(N,1);
        DQA1F12100h_P = zeros(N,1);
    elseif toggle_4>=1 && toggle_4<=5    % avg annual precip
        DQA2Ensemble2100h_P = DeltaA2Ensemble2100_OW(:,1) * BetaHomog(S+KCG+1);
        DQA1F12100h_P = DeltaA1F12100_OW(:,1) * BetaHomog(S+KCG+1);
    else                            % summer and winter precip
        DQA2Ensemble2100h_P1 = DeltaA2Ensemble2100_OW3(:,1) * BetaHomog(S+KCG+1);
        DQA2Ensemble2100h_P2 = DeltaA2Ensemble2100_OW3(:,2) * BetaHomog(S+KCG+2);
        DQA2Ensemble2100h_P = DQA2Ensemble2100h_P1 + DQA2Ensemble2100h_P2;
        DQA1F12100h_P1 = DeltaA1F12100_OW3(:,1) * BetaHomog(S+KCG+1);
        DQA1F12100h_P2 = DeltaA1F12100_OW3(:,2) * BetaHomog(S+KCG+2);
        DQA1F12100h_P = DQA1F12100h_P1 + DQA1F12100h_P2;    end
    AggWelfareA2Ensemble2100h_P = DQA2Ensemble2100h_P' * Income;          % in billions
    WelfarePctA2Ensemble2100h_P = AggWelfareA2Ensemble2100h_P / sum(Income);
    AggWelfareA2Ensemble2100h_P = WelfarePctA2Ensemble2100h_P * USPersIncome2000;
    AggWelfareA1F12100h_P = DQA1F12100h_P' * Income;          % in billions
    WelfarePctA1F12100h_P = AggWelfareA1F12100h_P / sum(Income);
    AggWelfareA1F12100h_P = WelfarePctA1F12100h_P * USPersIncome2000;
    
    % Welfare impact of humid change
    if toggle_4<2        % humid not in model
        DQA2Ensemble2100h_H = zeros(N,1);
        DQA1F12100h_H = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=3    % average humid
        DQA2Ensemble2100h_H = DeltaA2Ensemble2100_OW(:,2) * BetaHomog(S+KCG+2);
        DQA1F12100h_H = DeltaA1F12100_OW(:,2) * BetaHomog(S+KCG+2);
    elseif toggle_4>=4 && toggle_4<=5    % summer humid (June-Aug)
        DQA2Ensemble2100h_H = DeltaA2Ensemble2100_OW2(:,2) * BetaHomog(S+KCG+2);
        DQA1F12100h_H = DeltaA1F12100_OW2(:,2) * BetaHomog(S+KCG+2);        
    else                            % summer and winter humid
        DQA2Ensemble2100h_H1 = DeltaA2Ensemble2100_OW3(:,3) * BetaHomog(S+KCG+3);
        DQA2Ensemble2100h_H2 = DeltaA2Ensemble2100_OW3(:,4) * BetaHomog(S+KCG+4);
        DQA2Ensemble2100h_H = DQA2Ensemble2100h_H1 + DQA2Ensemble2100h_H2;
        DQA1F12100h_H1 = DeltaA1F12100_OW3(:,3) * BetaHomog(S+KCG+3);
        DQA1F12100h_H2 = DeltaA1F12100_OW3(:,4) * BetaHomog(S+KCG+4);
        DQA1F12100h_H = DQA1F12100h_H1 + DQA1F12100h_H2;        
    end
    AggWelfareA2Ensemble2100h_H = DQA2Ensemble2100h_H' * Income;          % in billions
    WelfarePctA2Ensemble2100h_H = AggWelfareA2Ensemble2100h_H / sum(Income);
    AggWelfareA2Ensemble2100h_H = WelfarePctA2Ensemble2100h_H * USPersIncome2000;
    AggWelfareA1F12100h_H = DQA1F12100h_H' * Income;          % in billions
    WelfarePctA1F12100h_H = AggWelfareA1F12100h_H / sum(Income);
    AggWelfareA1F12100h_H = WelfarePctA1F12100h_H * USPersIncome2000;
    
    % Welfare impact of sunshine change
    % Be careful here. Control var in sunshine fraction but model change is
    % cloud fraction. So need to take negative.
    if toggle_4<2           % sunshine not in model
        DQA2Ensemble2100h_S = zeros(N,1);
        DQA1F12100h_S = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=5              % average annual sunshine
        DQA2Ensemble2100h_S = -DeltaA2Ensemble2100_OW(:,3) * BetaHomog(S+KCG+3);
        DQA1F12100h_S = -DeltaA1F12100_OW(:,3) * BetaHomog(S+KCG+3);
    else                            % summer and winter sunshine
        DQA2Ensemble2100h_S1 = -DeltaA2Ensemble2100_OW3(:,5) * BetaHomog(S+KCG+5);
        DQA2Ensemble2100h_S2 = -DeltaA2Ensemble2100_OW3(:,6) * BetaHomog(S+KCG+6);
        DQA2Ensemble2100h_S = DQA2Ensemble2100h_S1 + DQA2Ensemble2100h_S2;
        DQA1F12100h_S1 = -DeltaA1F12100_OW3(:,5) * BetaHomog(S+KCG+5);
        DQA1F12100h_S2 = -DeltaA1F12100_OW3(:,6) * BetaHomog(S+KCG+6);
        DQA1F12100h_S = DQA1F12100h_S1 + DQA1F12100h_S2;
    end
    AggWelfareA2Ensemble2100h_S = DQA2Ensemble2100h_S' * Income;          % in billions
    WelfarePctA2Ensemble2100h_S = AggWelfareA2Ensemble2100h_S / sum(Income);
    AggWelfareA2Ensemble2100h_S = WelfarePctA2Ensemble2100h_S * USPersIncome2000;
    AggWelfareA1F12100h_S = DQA1F12100h_S' * Income;          % in billions
    WelfarePctA1F12100h_S = AggWelfareA1F12100h_S / sum(Income);
    AggWelfareA1F12100h_S = WelfarePctA1F12100h_S * USPersIncome2000;    
    
    % Put together overall welfare impact and set up output
    DQA2Ensemble2100h = DQA2Ensemble2100h_T + DQA2Ensemble2100h_P + DQA2Ensemble2100h_H + DQA2Ensemble2100h_S;
    AggWelfareA2Ensemble2100h = DQA2Ensemble2100h' * Income;          % in billions
    WelfarePctA2Ensemble2100h = AggWelfareA2Ensemble2100h / sum(Income);
    AggWelfareA2Ensemble2100h = WelfarePctA2Ensemble2100h * USPersIncome2000;   
    DQA1F12100h = DQA1F12100h_T + DQA1F12100h_P + DQA1F12100h_H + DQA1F12100h_S;
    AggWelfareA1F12100h = DQA1F12100h' * Income;          % in billions
    WelfarePctA1F12100h = AggWelfareA1F12100h / sum(Income);
    AggWelfareA1F12100h = WelfarePctA1F12100h * USPersIncome2000;  
    PctLosingA2h = sum(PopStored.*(DQA2Ensemble2100h<0))/sum(PopStored);        % pct of population losing
    PctLosingA1F1h = sum(PopStored.*(DQA1F12100h<0))/sum(PopStored);        % pct of population losing
    
    
    % FINAL PIECE: MIGRATION
    ES = 8;      % Housing supply elasticity
    WelfarePctA2Ensemble2100hPop = DQA2Ensemble2100h' * PopStored / sum(PopStored); % pop-weighted
    DN = ES * (DQA2Ensemble2100h - WelfarePctA2Ensemble2100hPop);  % change in log population
    IncomeNew = exp(log(Income) + DN);
    WelfarePctA2Ensemble2100h_Migration = DQA2Ensemble2100h' * IncomeNew / sum(IncomeNew);
    PctPopChg = exp(DN) - 1;    % percent change in population
    
    % Store point estimates and analytic std errors on WTP
    BinValueh40_365 = BinValueh(Bin40)*365;
    SDBinValueh40_365 = SDBinValueh(Bin40)*365;
    BinValueh80_365 = BinValueh(Bin80)*365;
    SDBinValueh80_365 = SDBinValueh(Bin80)*365;
    WelfarePctA2Ensemble2100h_TOtherWeatherComponent = -(WelfarePctA2Ensemble2100h_TW + WelfarePctA2Ensemble2100h_TS - WelfarePctA2Ensemble2100h);
    WelfarePctA1F12100h_TOtherWeatherComponent       = -(WelfarePctA1F12100h_TW       + WelfarePctA1F12100h_TS       - WelfarePctA1F12100h);
    ChangeWelfareinBillions2008A2Ensemble2100h = WelfarePctA2Ensemble2100h*USPersIncome2000;
    ChangeWelfareinBillions2008A1F12100h = WelfarePctA1F12100h*USPersIncome2000;
    
    HomogOutput = [BinValueh40_365 SDBinValueh40_365 BinValueh80_365 SDBinValueh80_365 WelfarePctA2Ensemble2100h WelfarePctA2Ensemble2100h_TW ];
    HomogOutput = [HomogOutput WelfarePctA2Ensemble2100h_TS WelfarePctA2Ensemble2100h_TOtherWeatherComponent ChangeWelfareinBillions2008A2Ensemble2100h]; 
    HomogOutput = [HomogOutput PctLosingA2h WelfarePctA1F12100h WelfarePctA1F12100h_TW WelfarePctA1F12100h_TS  WelfarePctA1F12100h_TOtherWeatherComponent ChangeWelfareinBillions2008A1F12100h PctLosingA1F1h]; 
    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % LOCAL LINEAR REGRESSION
    % LOOP OVER EACH OBSERVATION. RUN OLS, PUTTING WEIGHT ON SIMILAR
    % OBSERVATIONS. GET OBSERVATION-SPECIFIC COEFFICIENTS
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Identities of characteristics on which to allow random parameters
    KIndex = 1:K;
    KeepCols = KIndex(WCols==1); DropCols = KIndex(WCols==0);
    KP = length(KeepCols);
    % Set up vector of results
    BETA = zeros(N,KP); SDBETA = zeros(N,KP); Z = zeros(N,KP);      % initialize
    VARCOVAR = zeros(KP,KP,N);   VARCOVARNOCLUST = zeros(KP,KP,N);   % covariance matrix for BETA
    % Standard deviations of covariates
    Sigma = std(X(:,KeepCols));
    
    % Take effects of linear parameters / coefficients out of the price
    % variable
    if FE==0 && DFE==0
        QR = Q - X(:,DropCols) * BetaHomog(DropCols);
    elseif FE==0 && DFE==1   
        QR = Q - [X(:,DropCols) DivFE] * [BetaHomog(DropCols); BetaHomog(K+1:KKK)];
    else
        QR = Q - [X(:,DropCols) StateFE] * [BetaHomog(DropCols); BetaHomog(K+1:KKK)];
    end

    % Population weight dep var and selected regressors
    XR = diag(Pop.^0.5) * X(:,KeepCols);
    QR = diag(Pop.^0.5) * QR;

    % Set up characteristics to be used in calculation of kernel weights
    if toggle_12<=5 || toggle_12>=50            % use X(:,KeepCols)
        XW = X; KeepColsW = KeepCols; KPW = KP; SigmaW = Sigma;
    else                        % use lower-degree spline
        WColsW = WCols;
        WColsW(1:S) = 0;
        KeepColsW = KIndex(WColsW==1);     % only gets characteristics other than temperature spline
        XW = [XW X(:,KeepColsW)];       % combine low-order spline bases with other covariates
        KPW = size(XW,2);
        SigmaW = std(XW);
        KeepColsW = 1:KPW;              % Set to include all columns of XW
    end
    
    % Now loop over each observation, doing a local regression at each
    for i = 1:N
        Z = zeros(N,KPW);     % initialize 
        % Need to build up kernel weighting matrix for all observations
        % if mod(i,500)==0; disp(i); end;       % show progress every 200 obs
        W = ones(N,1);     % initialize weights
        for k = 2:KPW            % ignore constant, so start at 2
            kk = KeepColsW(k);
            Z(:,k) = (XW(:,kk) - XW(i,kk)) / SigmaW(k);  % distance from current obs
        end
        Z = Z / h;
        Z = normpdf(Z,0,1);
        for k = 1:KPW
            W = W .* Z(:,k);
        end
        W = (W / h).^0.5;     % weights (must take square root given back-division below)

        % Now do weighted least squares    
        QHat = diag(W) * QR;
        XHat = diag(W) * XR;
        Qinv = inv(XHat' * XHat);
        
        BETA(i,:) = [Qinv * XHat' * QHat]';    % local coefs at location i
        
        if toggle_11==1
            % Get variance of BETA(i)
            Resids = QHat - XHat * BETA(i,:)';
            XXe = XHat .* repmat(Resids,1,KP);
            VARCOVARNOCLUST(:,:,i) = N / (N-KKK) * (Qinv * XXe' * XXe * Qinv);
            % get clustered variance
            EE = zeros(N,N); 
            % EE = sparse(EE);
            for r = 1:N         % loop over rows of EE
                target = BS(r);     % msa of row r
                Ind = find(BS==target);     % all columns with msa equal to target
                EE(r,Ind) = 1;
                EE(r,:) = Resids(r) * Resids' .* EE(r,:);   % only capture residuals common to msa
            end
            VARCOVAR(:,:,i) = N / (N - KKK) * (Qinv * XHat' * EE * XHat * Qinv); % clustered variance
            clear EE
            SDBETA(i,:) = (diag(VARCOVAR(:,:,i)))'.^0.5;
        end
    end
	  
    % VARCOVAR = VARCOVARNOCLUST;  % UNCOMMENT TO USE NON-CLUSTERED VARIANCE
    
    
    % RESTORE FULL SET OF BETA PARAMETERS (USING HOMOG PREF ESTIMATES)
    BETAALL = repmat(BetaHomog',N,1);       % N by KKK
    BETAALL(:,KeepCols) = BETA;             % Insert heterog prefs
    BETA = BETAALL;
    clear BETAALL
    
    

    % Get value of an extra day at 65 degrees. Will need to subtract this off
    % of all estimates
    Val65 = BETA(:,1:S) * Basis(Bin65Ind,:)';

    % Get value of an extra day at each temperature in each location (rel to 65
    % degrees)
    BinValue = BETA(:,1:S) * Basis' - repmat(Val65,1,NbinInd);   % impact of one more day in each bin

    if toggle_11==1
        % Std error of marginal impact
        BasisDiff = Basis - repmat(Basis(Bin65Ind,:),NbinInd,1);
        VarBinValue = zeros(N,NbinInd);       % initialize
        for i = 1:N
            for b = 1:NbinInd
                VarBinValue(i,b) = BasisDiff(b,:) * VARCOVAR(1:S,1:S,i) * BasisDiff(b,:)';
            end
        end
        SDBinValue = VarBinValue.^0.5;
        TopInt = BinValue + 1.96 * SDBinValue;   % top of 95% c.i.
        BotInt = BinValue - 1.96 * SDBinValue;   % bottom of 95% c.i.
    end
        
    % If using splines, extrapolate MWTP at extremes
    if toggle_12>3 && toggle_12<50
        % Initialize valuations; BinValueh within present-day temp range, zero outside
        BinValueInd = BinValue;     % defined over present-day temp range
        BinValue = zeros(N,Nbin);   % defined over present & future temps (wider) 
        BinValue(:,Inds) = BinValueInd;
        if toggle_11==1      % do conf interval too
            TopIntInd = TopInt; BotIntInd = BotInt;
            TopInt = zeros(N,Nbin); BotInt = zeros(N,Nbin);
            TopInt(:,Inds) = TopIntInd; BotInt(:,Inds) = BotIntInd;
        end
        % Loop over Pumas, flattening valuations outside of the active
        % spline range
        for i = 1:N
            CenterS = CenterSh; LowS = LowSh; HighS = HighSh; MinS = MinSh; MaxS = MaxSh;
            ColdVal = BinValue(i,MinS); HotVal = BinValue(i,MaxS);    % valuations at bottom and top of temp dist
            BinValue(i,LowS) = ColdVal; BinValue(i,HighS) = HotVal;
            if toggle_11==1      % do conf interval too
                TopInt(i,LowS) = TopInt(i,MinS); TopInt(i,HighS) = TopInt(i,MaxS);
                BotInt(i,LowS) = BotInt(i,MinS); BotInt(i,HighS) = BotInt(i,MaxS);
            end    
        end
    end
        
    
    % TEST OF H0: MODEL IS HDD/CDD, NOT 5 PART LINEAR SPLINE
    if toggle_11==1
        if toggle_12==1
            TESTSTAT3 = zeros(N,1); TESTPVAL3 = zeros(N,1);
            for i = 1:N
                J = 2;
                R = [0 1 -1 0 0; 0 0 0 1 -1];
                TESTSTAT3(i) = (R*BETA(i,1:5)')' * inv(R * VARCOVAR(1:5,1:5,i) * R') * R*BETA(i,1:5)' / J;
                DF = N - K;
                TESTPVAL3(i) = 1 - fcdf(TESTSTAT3(i),J,DF);
            end
            TESTPVAL3_Med = median(TESTPVAL3);
            TESTPVAL3_Rej = length(find(TESTPVAL3<0.05));    
        end
        % TEST OF H0: MODEL IS FLAT AT EXTREMES, NOT 5 PART LINEAR SPLINE
        if toggle_12==1
            TESTSTAT0 = zeros(N,1); TESTPVAL0 = zeros(N,1);
            for i = 1:N
                J = 2;
                R = [0 1 0 0 0; 0 0 0 0 1];
                TESTSTAT0(i) = (R*BETA(i,1:5)')' * inv(R * VARCOVAR(1:5,1:5,i) * R') * R*BETA(i,1:5)' / J;
                DF = N - K;
                TESTPVAL0(i) = 1 - fcdf(TESTSTAT0(i),J,DF);
            end
            TESTPVAL0_Med = median(TESTPVAL0);
            TESTPVAL0_Rej = length(find(TESTPVAL0<0.05));    
        end
    end
    
   
    % Plot impact of an extra day in a selected bin
    if toggle_11==1
        % First get CDD, HDD variables
        Helper = zeros(Nbin,1);
        for i = Bin65:Nbin
            Helper(i) = i - Bin65;     % line on CDD side
        end        
        CDD = Bin * Helper * 0.9;
        Helper = zeros(Nbin,1);
        for i = 1:Bin65
            Helper(i) = Bin65 - i;     % line on HDD side
        end 
        HDD = Bin * Helper * 0.9;
        BinNum = floor((40-Tmin+0.449)/0.9)+1;  % select bin (enter desired temperature in F)
        %Char = Bin(:,BinNum);          % # of days in bin
        %Char = XCOW(:,2);               % to plot against summer humid instead
        Char = HDD;
        %Char = CDD;
        Coef = BinValue(:,BinNum) * 365;     % impact of one more day in bin
        if toggle_12>3 && toggle_12<50
            BinNumSD = floor((40-TminInd+0.449)/0.9)+1;
        else
            BinNumSD = floor((40-Tmin+0.449)/0.9)+1;
        end
        SDCoef = SDBinValue(:,BinNumSD) * 365;   % s.d. of impact
        %Coef = BETA(:,2);
        clf
        plot(Char,Coef,'x'); grid on
        xlabel('Heating degree days'); ylabel('MWTP for a day at 40F')
        %xlabel('Cooling degree days'); ylabel('MWTP for a day at 80F')
        axis([0 12000 -0.3 0])
        %axis([0 4500 -0.3 0])   
        % axis([0 4500 -0.8 0.2])    % For set ID
        plot(Char,SDCoef,'x'); grid on
        xlabel('Heating degree days'); ylabel('Std Dev of MWTP for a day at 40F')
        %xlabel('Cooling degree days'); ylabel('Std Dev of MWTP for a day at 80F')
        axis([0 12000 0 0.14])
        %axis([0 4500 0 0.14])         
        
        
        % Plot marginal impact function at a particular puma
        PUMA = 1703510;
        Ind = find(Puma==PUMA);       % location of the Puma code
        func = BinValue(Ind,:);             % marginal impact at all temperatures
        funchigh = TopInt(Ind,:); funclow = BotInt(Ind,:);
        tempminplot = Temps(Bin0); tempmaxplot = Temps(Bin110);
        temp = [tempminplot:0.9:tempmaxplot];  % vector of temperatures to plot
        Daycount = Bin(Ind,Bin0:Bin110);        % # of days in each bin
        DaycountA2Ensemble2100 = BinA2Ensemble2100(Ind,Bin0:Bin110);
        DaycountA1F12100 = BinA1F12100(Ind,Bin0:Bin110);
        clf
        [AX,H1,H2] = plotyy(temp,func(Bin0:Bin110),temp,Daycount);
        set(AX(1),'YLim',[-0.002 0.0005]); set(AX(2),'YLim',[0 20]);
        y1inc = (0.002 + 0.0005) / 5; y2inc = (20 - 0) / 5;
        set(AX(1),'YTick',[-0.002:y1inc:0.0005]); set(AX(2),'YTick',[0:y2inc:20]);
        grid on
        set(H1,'Color','b'); set(H1,'LineWidth',3);      % preference estimate
        set(H2,'LineStyle','--');  set(H2,'LineWidth',2); set(H2,'Color','k');  % current temps
        axes(AX(1)); ylabel('Marginal Value of an Extra Day at Each Temperature');
        axes(AX(2)); ylabel('Number of Days at Each Temperature');
        % add 2100 temps, top of c.i., and bottom of c.i.
        H3 = line(temp,DaycountA2Ensemble2100,'LineStyle','-.','LineWidth',2,'Color','r','Parent',AX(2));
        H4 = line(temp,funchigh(Bin0:Bin110),'LineStyle','--','LineWidth',2,'Color','b','Parent',AX(1));
        H5 = line(temp,funclow(Bin0:Bin110),'LineStyle','--','LineWidth',2,'Color','b','Parent',AX(1));
        set(AX(1),'XLim',[0 110]); set(AX(1),'XTick',[0:10:110]);
        set(AX(2),'XLim',[0 110]); set(AX(2),'XTick',[0:10:110]);
    end
    
    % Get weighted average MWTP across all PUMAs
    BinValueAvg = Income' * BinValue;
    BinValueAvg = BinValueAvg / sum(Income);
    BinValueAvg = BinValueAvg';
    
    
    
    % CONCAVITY CALCULATIONS
    SigmaX = std(XAll);    % std dev of each X variable (unweighted, including FE). For scaling   
    % Save the results of concavity of utility check in variable A
    Aij=zeros(N,N);     % Concavity needed to stop person at i from preferring j
    % Identify variables over which prefs are concave: all those that allow
    % random coefs (don't include constant term in linear spline models)
    if toggle_12<=3 || toggle_12>=50
        KeepColsC = KeepCols(2:KP);
    else
        KeepColsC = KeepCols(1:KP);     % cubic spline
    end
    % Get vector of unobserved attribute at each location
    Xi = zeros(N,1);
    for i = 1:N
        Xi(i) = Q(i) - BETA(i,:) * XAll(i,:)';
    end
    % Loop over i to get Aij. Numerator is linear utility difference
    % between j and i. Denominator is quadratic characteristic difference
    for i = 1:N
        numerator = Q(i) - Q + Xi - Xi(i);
        numerator = numerator' + BETA(i,:) * XAll' - BETA(i,:) * XAll(i,:)';    % 1 X N vector
        sum_quadratic_term = sum(((ones(N,1)*XAll(i,KeepColsC)-XAll(:,KeepColsC)) ./ (ones(N,1)*SigmaX(KeepColsC))).^2, 2);
        sum_quadratic_term = sum_quadratic_term.^0.5;
        Aij(i,:) = numerator ./sum_quadratic_term';  % concavity needed to stop i from preferring j
        Aij(i,i) = 0; % Since denominator is 0 on the diagonals.
    end
    MaxConcav = max([Aij zeros(N,1)],[],2); % N by 1. Concavity needed at each i to stop i from preferring all other locations    
    
    
    
    BETAPOINT = BETA;       % Store point-identified preferences before doing set identification
    BinValuePoint = BinValue;
    if toggle_10==1
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % SET IDENTIFICATION OF PREFERENCES
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Follows the Bajari and Benkard NBER WP version. Gibbs sampler to get
        % distribution of possible preference parameters at each location

        % Start by obtaining the parameter space. Increase the size of the space
        % in each dimension relative to what was found with point
        % identification
        ThetaRange = max(BETA(:,1:K)) - min(BETA(:,1:K));
        ThetaMin = min(BETA(:,1:K)) - ThetaRange/2;
        ThetaMax = max(BETA(:,1:K)) + ThetaRange/2;
        % Set number of draws
        D = 1500;
        DMin = 501;     % first draw to keep
        % Set matrix of parameter draws
        BETASET = zeros(N,K,D-DMin+1);

        Counter = 1:N;      % handy index of all observations
        
        % Loop over all locations
        for i = 1:N
            % Initalize results matrix B
            B = zeros(D,KKK);
            B(:,1) = BETA(i,1);             % Initialize constant for all draws d
            B(:,K+1:KKK) = ones(D,1)*BETA(i,K+1:KKK); % Initialize fixed effects for all draws d
            BMin = zeros(D,K); BMax = zeros(D,K);               % initialize
            for d = 1:D             % draw D Beta's at current location i
                % Initial seed: use point identified estimate
                if d==1; 
                    B(d,:) = BETA(i,:); % stored set of coefficients
                    BC = BETA(i,:);     % active set of coefficients
                    % See if guess is consistent with utility maximization
                    % 'test' variable must be greater than or equal to zero
                    sum_quadratic_term = sum(((ones(N,1)*XAll(i,KeepColsC)-XAll(:,KeepColsC)) ./ (ones(N,1)*SigmaX(KeepColsC))).^2, 2);
                    sum_quadratic_term = sum_quadratic_term.^0.5;
                    test = BC*XAll(i,:)' - XAll*BC' + MaxConcav(i)*sum_quadratic_term;
                    test = test + Q - Q(i) + Xi(i) - Xi;
                    Test = min(test(Counter~=i));
                else
                end
                % BEGIN GIBBS SAMPLING ALGORITHM
                % LOOP THROUGH PARAMETERS ONE BY ONE (EXCEPT CONSTANT)
                % FIRST OBTAIN DISTRIBUTION, CONDITIONAL ON OTHER PARAMETER VALUES
                % AND CHOOSING LOCATION i
                % SECOND DRAW FROM THAT DISTRIBUTION
                for k = 2:K
                    if WCols(k)==0      % coefficient k is constant at all locations
                        BC(k) = BETA(i,k);  % use point identified estimate
                        B(d,k) = BC(k);
                    else
                        % Take other coefs as given when forming inequality to get
                        % bounds of distribution
                        RHS = XAll(:,1:k-1)*BC(1:k-1)' - BC(1:k-1)*XAll(i,1:k-1)';
                        if k<K
                            RHS = RHS + XAll(:,k+1:K)*BC(k+1:K)' - BC(k+1:K)*XAll(i,k+1:K)';
                        else
                        end
                        % Add fixed effects if a fixed effect model
                        if FE==1 || DFE==1
                            RHS = RHS + XAll(:,K+1:KKK)*BC(K+1:KKK)' - BC(1,K+1:KKK)*XAll(i,K+1:KKK)';
                        else
                        end
                        % Add in quadratic term, unobservable, and prices
                        RHS = RHS + Xi - Xi(i) + Q(i) - Q - MaxConcav(i)*sum_quadratic_term;
                        % Each j gives either an upper or lower bound on Beta_ik. Split
                        % vector by sign of difference in characteristic k
                        GL = XAll(i,k) - XAll(:,k);   % sign determines whether we look for upper or lower bound
                        IndLB = find(GL>0);     % j's that yield a lower bound
                        IndUB = find(GL<0);     % j's that yield an upper bound
                        if isempty(IndLB)       % characteristic k at location i is extremely low
                            BLB = ThetaMin(k);  % lower bound of parameter space
                        else
                            BLB = RHS(IndLB) ./ (XAll(i,k)-XAll(IndLB,k)); 
                        end
                        if isempty(IndUB)       % characteristic k at location i is extremely high
                            BUB = ThetaMax(k);  % upper bound of parameter space
                        else
                            BUB = RHS(IndUB) ./ (XAll(i,k)-XAll(IndUB,k)); 
                        end
                        BMin(d,k) = max([BLB; ThetaMin(k)]); % satisfy all inequalities and stay in parameter space
                        BMax(d,k) = min([BUB; ThetaMax(k)]); % satisfy all inequalities and stay in parameter space
                        % Now draw a parameter from this distribution and store it
                        BC(k) = BMin(d,k) + rand() * (BMax(d,k)-BMin(d,k));
                        B(d,k) = BC(k);
                    end
                end
            end
            BETASET(i,:,:) = B(DMin:D,1:K)';        % store results
            if mod(i,100)==0; disp(i); else end     % show progress
        end


        % Overall comparison
        mean(mean(squeeze(BETASET(:,3,:))'));
        mean(BETA(:,3));

        % GET MEAN PREFERENCES AT EACH LOCATION; USE THE MEANS TO DO WELFARE
        % RESULTS ARE VALID WITH WELFARE FUNCTIONS THAT ARE LINEAR IN
        % PARAMETERS AND IGNORING MIGRATION
        if FE==1 || DFE==1
            BETA = [mean(BETASET,3) BETAPOINT(:,K+1:KKK)]; % mean of set-identified preferences (distribution is point-identified)
        else
            BETA = mean(BETASET,3);  % mean of set-identified preferences (distribution is point-identified)
        end    
        BETASETAVG = BETA;      % store mean of set-identified preferences
        
        % Get value of an extra day at 65 degrees. Will need to subtract this off
        % of all estimates
        Val65 = BETA(:,1:S) * Basis(Bin65Ind,:)';
        
        % Get value of an extra day at each temperature in each location (rel to 65
        % degrees)
        BinValue = BETA(:,1:S) * Basis' - repmat(Val65,1,NbinInd);   % impact of one more day in each bin    
        
        % Get values at 40, 65, and 80F at all locations and draws
        Val65SET = zeros(N,D-DMin+1);
        Val40SET = zeros(N,D-DMin+1);
        Val80SET = zeros(N,D-DMin+1);
        for i = 1:N
            Val65SET(i,:) = (squeeze(BETASET(i,1:S,:))' * Basis(Bin65Ind,:)')';
            Val40SET(i,:) = (squeeze(BETASET(i,1:S,:))' * Basis(Bin40Ind,:)')' - Val65SET(i,:);
            Val80SET(i,:) = (squeeze(BETASET(i,1:S,:))' * Basis(Bin80Ind,:)')' - Val65SET(i,:);
        end
        
        
        % SET UP EXPORT FILE TO STATA FOR KERNEL DENSITY PLOT
        % First get CDD, HDD variables
        Helper = zeros(Nbin,1);
        for i = Bin65:Nbin
            Helper(i) = i - Bin65;     % line on CDD side
        end        
        CDD = Bin * Helper * 0.9;
        Helper = zeros(Nbin,1);
        for i = 1:Bin65
            Helper(i) = Bin65 - i;     % line on HDD side
        end 
        HDD = Bin * Helper * 0.9;
        ExportCold = [Puma Income HDD BinValuePoint(:,Bin40) BinValue(:,Bin40Ind) Val40SET];
        ExportHeat = [Puma Income CDD BinValuePoint(:,Bin80) BinValue(:,Bin80Ind) Val80SET];
        dlmwrite('SetID\ExportCold_Dem2FE0_TempWtOnly_S7_99_h2_2x.csv',ExportCold,'precision',8);
        dlmwrite('SetID\ExportHeat_Dem2FE0_TempWtOnly_S7_99_h2_2x.csv',ExportHeat,'precision',8);
        save SetID\BETASET_Dem2FE0_TempWtOnly_S7_99_h2_2x.mat BETASET
        
        blah = BinValue;
        
        % If using splines, extrapolate MWTP at extremes
        if toggle_12>3 && toggle_12<50
            % Initialize valuations; BinValueh within present-day temp range, zero outside
            BinValueInd = BinValue;     % defined over present-day temp range
            BinValue = zeros(N,Nbin);   % defined over present & future temps (wider) 
            BinValue(:,Inds) = BinValueInd;
            % Loop over Pumas, flattening valuations outside of the active
            % spline range
            for i = 1:N
                CenterS = CenterSh; LowS = LowSh; HighS = HighSh; MinS = MinSh; MaxS = MaxSh;
                ColdVal = BinValue(i,MinS); HotVal = BinValue(i,MaxS);    % valuations at bottom and top of temp dist
                BinValue(i,LowS) = ColdVal; BinValue(i,HighS) = HotVal;  
            end
        end
        % Get weighted average MWTP across all PUMAs
        BinValueAvgPoint = BinValueAvg;     % store point identified MWTPs
        BinValueAvg = Income' * BinValue;
        BinValueAvg = BinValueAvg / sum(Income);
        BinValueAvg = BinValueAvg';  
    else                        % Set identification toggle_10 is off
    end

    
    
    
    
    %%%%%%%%%%%%
    % HETEROGENEOUS PREF WELFARE CALCULATIONS
    % Start with welfare impact of temperature change; linear terms
    DQA2Ensemble2100_T = zeros(N,1);     % initialize
    DQA1F12100_T = zeros(N,1);     % initialize
    for i = 1:N
        DQA2Ensemble2100_T(i) = BinValue(i,:) * DeltaA2Ensemble2100(i,:)';     % delta welfare as pct of income
        DQA1F12100_T(i) = BinValue(i,:) * DeltaA1F12100(i,:)';     % delta welfare as pct of income
    end
    
    % Additional Concavity Term. All values are a percent of income
    if toggle_12>3 && toggle_12<50      % Flatten out extremes of spline and create new X's
        BasisF = zeros(Nbin,S);       % initalize full basis
        BasisF(Inds,:) = Basis;       % basis functions over present-day temperatures
        BasisF(LowS,:) = repmat(BasisF(MinS,:),length(LowS),1); % flatten basis at cold temps
        BasisF(HighS,:) = repmat(BasisF(MaxS,:),length(HighS),1); % flatten basis at cold temps
        XAllF(:,1:S) = Bin * BasisF;  % covariates corrected for flattened spline
        SigmaXF = std(XAllF);    % std dev of each X variable. For scaling   
        NonZeroXF = find(SigmaXF>0);        % covariates that are non-zero after flattening
        XAllF = XAllF(:,NonZeroXF); SigmaXF = SigmaXF(NonZeroXF);
        NonZeroSF = find(NonZeroXF<=S);       % elements of NonZeroX that pertain to splines
        BasisF = BasisF(:,NonZeroXF(NonZeroSF));  % non-degenerate basis functions
        SF = size(BasisF,2);            % number of remaining basis functions
    else
        SF = S; BasisF = Basis; XAllF = XAll; SigmaXF = SigmaX;
    end
    % Quadratic welfare losses
    % Include constant only if cubic spline model
    if toggle_12<=3 || toggle_12>=50; s=2; else; s=1; end;
    DQQA2Ensemble2100 = - sum(((BinA2Ensemble2100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
    DQQA1F12100 = - sum(((BinA1F12100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;

    % Get change in marginal WTP given concavity and A2 ensemble warming
    BinValueQ = zeros(N,Nbin);      % initialize
    diff65 = (BinA2Ensemble2100-Bin)*BasisF.*repmat(BasisF(Bin65,:),N,1);   % N by SF. Change at 65F
    for b = 1:Nbin
        diff = (BinA2Ensemble2100-Bin)*BasisF.*repmat(BasisF(b,:),N,1);     % N by SF. Change at bin b
        BinValueQ(:,b) = -2*MaxConcav .* sum(((diff-diff65)./repmat(SigmaXF(1:SF).^2,N,1)),2);
    end
    clear diff diff65
        
    
    % Aggregate welfare in billions. Linear and quadratic terms
    AggWelfareA2Ensemble2100_T = DQA2Ensemble2100_T' * Income;          % in billions
    AggWelfareA1F12100_T = DQA1F12100_T' * Income;          % in billions
    AggWelfareQA2Ensemble2100 = DQQA2Ensemble2100' * Income;          % in billions
    AggWelfareQA1F12100 = DQQA1F12100' * Income;          % in billions
   
    
    % as a percent of total income. Linear and quadratic terms
    WelfarePctA2Ensemble2100_T = AggWelfareA2Ensemble2100_T / sum(Income);
    WelfarePctA1F12100_T = AggWelfareA1F12100_T / sum(Income);
    WelfarePctQA2Ensemble2100 = AggWelfareQA2Ensemble2100 / sum(Income);
    WelfarePctQA1F12100 = AggWelfareQA1F12100 / sum(Income);

    
    % Get welfare impact in $ using total US personal income
    AggWelfareA2Ensemble2100_T = WelfarePctA2Ensemble2100_T * USPersIncome2000;
    AggWelfareA1F12100_T = WelfarePctA1F12100_T * USPersIncome2000;
    AggWelfareQA2Ensemble2100 = WelfarePctQA2Ensemble2100 * USPersIncome2000;
    AggWelfareQA1F12100 = WelfarePctQA1F12100 * USPersIncome2000;

    
    
    
    
    
    % Welfare impact of precip change
    if toggle_4<1       % precipitation not in model
        DQA2Ensemble2100_P = zeros(N,1);
        DQA1F12100_P = zeros(N,1);
    elseif toggle_4>=1 && toggle_4<=5    % avg annual precip
        DQA2Ensemble2100_P = DeltaA2Ensemble2100_OW(:,1) .* BETA(:,S+KCG+1);
        DQA1F12100_P = DeltaA1F12100_OW(:,1) .* BETA(:,S+KCG+1);
    else                            % summer and winter precip
        DQA2Ensemble2100_P1 = DeltaA2Ensemble2100_OW3(:,1) .* BETA(:,S+KCG+1);
        DQA2Ensemble2100_P2 = DeltaA2Ensemble2100_OW3(:,2) .* BETA(:,S+KCG+2);
        DQA2Ensemble2100_P = DQA2Ensemble2100h_P1 + DQA2Ensemble2100h_P2;
        DQA1F12100_P1 = DeltaA1F12100_OW3(:,1) .* BETA(:,S+KCG+1);
        DQA1F12100_P2 = DeltaA1F12100_OW3(:,2) .* BETA(:,S+KCG+2);
        DQA1F12100_P = DQA1F12100h_P1 + DQA1F12100h_P2;
    end
    AggWelfareA2Ensemble2100_P = DQA2Ensemble2100_P' * Income;          % in billions
    WelfarePctA2Ensemble2100_P = AggWelfareA2Ensemble2100_P / sum(Income);
    AggWelfareA2Ensemble2100_P = WelfarePctA2Ensemble2100_P * USPersIncome2000;
    AggWelfareA1F12100_P = DQA1F12100_P' * Income;          % in billions
    WelfarePctA1F12100_P = AggWelfareA1F12100_P / sum(Income);
    AggWelfareA1F12100_P = WelfarePctA1F12100_P * USPersIncome2000;
    
    % Welfare impact of summer humid change
    if toggle_4<2        % humid not in model
        DQA2Ensemble2100_H = zeros(N,1);
        DQA1F12100_H = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=3    % average annual humid
        DQA2Ensemble2100_H = DeltaA2Ensemble2100_OW(:,2) .* BETA(:,S+KCG+2);
        DQA1F12100_H = DeltaA1F12100_OW(:,2) .* BETA(:,S+KCG+2);
    elseif toggle_4>=4 && toggle_4<=5    % summer humid
        DQA2Ensemble2100_H = DeltaA2Ensemble2100_OW2(:,2) .* BETA(:,S+KCG+2);
        DQA1F12100_H = DeltaA1F12100_OW2(:,2) .* BETA(:,S+KCG+2);
    else                            % summer and winter humid
        DQA2Ensemble2100_H1 = DeltaA2Ensemble2100_OW3(:,3) .* BETA(:,S+KCG+3);
        DQA2Ensemble2100_H2 = DeltaA2Ensemble2100_OW3(:,4) .* BETA(:,S+KCG+4);
        DQA2Ensemble2100_H = DQA2Ensemble2100h_H1 + DQA2Ensemble2100h_H2;
        DQA1F12100_H1 = DeltaA1F12100_OW3(:,3) .* BETA(:,S+KCG+3);
        DQA1F12100_H2 = DeltaA1F12100_OW3(:,4) .* BETA(:,S+KCG+4);
        DQA1F12100_H = DQA1F12100h_H1 + DQA1F12100h_H2;
    end
    AggWelfareA2Ensemble2100_H = DQA2Ensemble2100_H' * Income;          % in billions
    WelfarePctA2Ensemble2100_H = AggWelfareA2Ensemble2100_H / sum(Income);
    AggWelfareA2Ensemble2100_H = WelfarePctA2Ensemble2100_H * USPersIncome2000;
    AggWelfareA1F12100_H = DQA1F12100_H' * Income;          % in billions
    WelfarePctA1F12100_H = AggWelfareA1F12100_H / sum(Income);
    AggWelfareA1F12100_H = WelfarePctA1F12100_H * USPersIncome2000;
    
    % Welfare impact of sunshine change
    % Be careful here. Control var in sunshine fraction but model change is
    % cloud fraction. So need to take negative.
    if toggle_4<2           % sunshine not in model
        DQA2Ensemble2100_S = zeros(N,1);
        DQA1F12100_S = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=5              % average annual sunshine
        DQA2Ensemble2100_S = -DeltaA2Ensemble2100_OW(:,3) .* BETA(:,S+KCG+3);
        DQA1F12100_S = -DeltaA1F12100_OW(:,3) .* BETA(:,S+KCG+3);
    else                            % summer and winter sunshine
        DQA2Ensemble2100_S1 = -DeltaA2Ensemble2100_OW3(:,5) .* BETA(:,S+KCG+5);
        DQA2Ensemble2100_S2 = -DeltaA2Ensemble2100_OW3(:,6) .* BETA(:,S+KCG+6);
        DQA2Ensemble2100_S = DQA2Ensemble2100h_S1 + DQA2Ensemble2100h_S2;
        DQA1F12100_S1 = -DeltaA1F12100_OW3(:,5) .* BETA(:,S+KCG+5);
        DQA1F12100_S2 = -DeltaA1F12100_OW3(:,6) .* BETA(:,S+KCG+6);
        DQA1F12100_S = DQA1F12100h_S1 + DQA1F12100h_S2;
    end
    AggWelfareA2Ensemble2100_S = DQA2Ensemble2100_S' * Income;          % in billions
    WelfarePctA2Ensemble2100_S = AggWelfareA2Ensemble2100_S / sum(Income);
    AggWelfareA2Ensemble2100_S = WelfarePctA2Ensemble2100_S * USPersIncome2000;
    AggWelfareA1F12100_S = DQA1F12100_S' * Income;          % in billions
    WelfarePctA1F12100_S = AggWelfareA1F12100_S / sum(Income);
    AggWelfareA1F12100_S = WelfarePctA1F12100_S * USPersIncome2000;
    
    % Put together overall welfare impact and set up output
    DQA2Ensemble2100 = DQA2Ensemble2100_T + DQA2Ensemble2100_P + DQA2Ensemble2100_H + DQA2Ensemble2100_S;
    AggWelfareA2Ensemble2100 = DQA2Ensemble2100' * Income;          % in billions
    WelfarePctA2Ensemble2100 = AggWelfareA2Ensemble2100 / sum(Income);
    AggWelfareA2Ensemble2100 = WelfarePctA2Ensemble2100 * USPersIncome2000;   
    DQA1F12100 = DQA1F12100_T + DQA1F12100_P + DQA1F12100_H + DQA1F12100_S;
    AggWelfareA1F12100 = DQA1F12100' * Income;          % in billions
    WelfarePctA1F12100 = AggWelfareA1F12100 / sum(Income);
    AggWelfareA1F12100 = WelfarePctA1F12100 * USPersIncome2000;  
    PctLosingA2 = sum(PopStored.*(DQA2Ensemble2100<0))/sum(PopStored);        % pct of population losing
    PctLosingA1F1 = sum(PopStored.*(DQA1F12100<0))/sum(PopStored);        % pct of population losing    
    
    % Store point estimates
    BinValueAvg_Bin40_365 = BinValueAvg(Bin40)*365; 
    BinValueAvg_Bin80_365 =  BinValueAvg(Bin80)*365;
    
    LinearOutput = [BinValueAvg_Bin40_365 BinValueAvg_Bin80_365 WelfarePctA2Ensemble2100 WelfarePctA1F12100];

    % Write to Output
    save ('PartONE_POINTEstimate_Q1_D2_FE0_S7_99.mat', 'LinearOutput', 'HomogOutput');
    save ('PartONE_ALL_Q1_D2_FE0_S7_99.mat');
    
    % Display output so far
    display('Homogenous preference point estimates:')
    HomogOutput
    display('Heterogeneous preference point estimates:')
    LinearOutput
    
    
    
    
    
 %{
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXTRA CODE TO PLOT TEMPERATURE DISTRIBUTIONS
PUMA = 2603200
%PUMA = 4804608;
%PUMA = 602201;
Ind = find(Puma==PUMA);       % location of the FIP code
Daycount = Bin(Ind,Bin0:Bin110);        % # of days in each bin
DaycountA2Ensemble2100 = BinA2Ensemble2100(Ind,Bin0:Bin110);
clf
plot(temp,Daycount,'k--','LineWidth',3)
hold all
plot(temp,DaycountA2Ensemble2100,'r-.','LineWidth',2)
axis([0 110 0 25])
ylabel('Number of days at each temperature'); xlabel('Temperature (Fahrenheit)')
grid
 %}
    

    
    
    
    
%***************************************************************%
% BOOTSTRAPPING FOR STANDARD ERRORS
%***************************************************************%

% Save original point estimates. Homogenous prefs then heterog
KKK = size(XXX,2);        % # of vars in regression
Qinvh = inv(XXX' * XXX);
BetaHomogOrig = BetaHomog;
XOrig = XAll;
ResidsHomogOrig = QQQ - XXX * BetaHomog;    % note: these are population-weighted

BETAOrig = BETA;
% Obtain residuals
ResidsOrig = zeros(N,1);
for i = 1:N
    ResidsOrig(i) = Q(i) - XOrig(i,:) * BETAOrig(i,:)';     % note: not population-weighted
end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % START OF BOOTSTRAPPING LOOP
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    toggle_10 = 0; toggle_11 = 0;      % force no set ID, no std errors or plots
    NRep = 200;       % # of bootstrap repetitions
    HOMOGOUTPUT = ones(NRep,12) * -999;
    LINEAROUTPUT = ones(NRep,4) * -999;
    NB = max(BS);       % total number of bootstrap groups (msas)
    
    
    for rep = 1:NRep
    
    % Now build up matrices of data to use in the regressions
    QQQb = zeros(N,1); Qb = zeros(N,1); % initialize dependent vars for homog and local linear regressions
    % Draw new dependent variables, for both homog and local linear regressions
    for i = 1:NB        % loop over number of msas / states
        Ind = find(BS==i);  % all observations for state i
        RWt = -1 + 2 * round(rand());       % 50/50 chance of -1 or 1      
        QQQb(Ind) = XXX(Ind,:) * BetaHomogOrig + RWt * ResidsHomogOrig(Ind);  % for homog regressions (pop weighted)
        % Now Qb for local linear regressions (not pop weighted; must do
        % this later in loop)
        for j = 1:length(Ind)
            Qb(Ind(j)) = XOrig(Ind(j),:) * BETAOrig(Ind(j),:)' + RWt * ResidsOrig(Ind(j));
        end
    end



    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % HOMOGENOUS PREFERENCE REGRESSION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    BetaHomog = Qinvh * XXX' * QQQb;
        
        
    % Get value of an extra day at 65 degrees. Will need to subtract this off
    % of all estimates
    Val65h = Basis(Bin65Ind,:) * BetaHomog(1:S);

    % Marginal impact
    BinValueh = Basis * BetaHomog(1:S) - repmat(Val65h,NbinInd,1);
 
        
    % If using splines, extrapolate MWTP at extremes
    if toggle_12>3 && toggle_12<50
        % Initialize valuations; BinValueh within present-day temp range, zero outside
        BinValuehInd = BinValueh;     % defined over present-day temp range
        BinValueh = zeros(Nbin,1);    % defined over present & future temps (wider)
        BinValueh(Inds') = BinValuehInd;
        if toggle_11==1      % do conf interval & std errors too
            TopInthInd = TopInth; BotInthInd = BotInth;
            TopInth = zeros(Nbin,1); BotInth = zeros(Nbin,1);
            TopInth(Inds') = TopInthInd; BotInth(Inds') = BotInthInd;
            SDBinValuehInd = SDBinValueh;
            SDBinValueh = zeros(Nbin,1);
            SDBinValueh(Inds') = SDBinValuehInd;
        end
        % Define range of temperatures over which splines are still part of MWTP
        Pctile = 1 - (1 - toggle_9) / 2;    % percent of distribution (income-weighted) to exclude on each side
        CenterSh = find(CumDays>(1-Pctile) & CumDays<Pctile);    % select center of distribution
        LowSh = find(CumDays<=(1-Pctile)); HighSh = find(CumDays>=Pctile);
        MinSh = min(CenterSh'); MaxSh = max(CenterSh'); 
        % Now flatten valuations outside of the active spline range
        ColdVal = BinValueh(MinSh); HotVal = BinValueh(MaxSh);    % valuations at bottom and top of temp dist
        BinValueh(LowSh) = ColdVal; BinValueh(HighSh) = HotVal;
        if toggle_11==1      % do conf interval too
            TopInth(LowSh) = TopInth(MinSh); TopInth(HighSh) = TopInth(MaxSh);
            BotInth(LowSh) = BotInth(MinSh); BotInth(HighSh) = BotInth(MaxSh);
        end
    end


    
    
    %%%%%%%%
    % HOMOGENOUS PREF WELFARE CALCULATIONS
    % Start with welfare impact of temperature change
    DQA2Ensemble2100h_T = DeltaA2Ensemble2100 *  BinValueh;   % vector of delta welfare as pct of income
    DQA1F12100h_T = DeltaA1F12100 *  BinValueh;   % vector of delta welfare as pct of income

    AggWelfareA2Ensemble2100h_T = DQA2Ensemble2100h_T' * Income;          % in billions
    AggWelfareA1F12100h_T = DQA1F12100h_T' * Income;          % in billions
    % as a percent of total income
    WelfarePctA2Ensemble2100h_T = AggWelfareA2Ensemble2100h_T / sum(Income);
    WelfarePctA1F12100h_T = AggWelfareA1F12100h_T / sum(Income);
    % Get welfare impact in $ using total US personal income
    AggWelfareA2Ensemble2100h_T = WelfarePctA2Ensemble2100h_T * USPersIncome2000;
    AggWelfareA1F12100h_T = WelfarePctA1F12100h_T * USPersIncome2000;
    
    % Break into winter gain vs. summer loss
    DQA2Ensemble2100h_TW = DeltaA2Ensemble2100(:,1:Bin65) *  BinValueh(1:Bin65);   % vector of delta welfare as pct of income
    DQA1F12100h_TW = DeltaA1F12100(:,1:Bin65) *  BinValueh(1:Bin65);   % vector of delta welfare as pct of income    
    DQA2Ensemble2100h_TS = DeltaA2Ensemble2100(:,Bin65+1:Nbin) *  BinValueh(Bin65+1:Nbin);   % vector of delta welfare as pct of income
    DQA1F12100h_TS = DeltaA1F12100(:,Bin65+1:Nbin) *  BinValueh(Bin65+1:Nbin);   % vector of delta welfare as pct of income    
    WelfarePctA2Ensemble2100h_TW = DQA2Ensemble2100h_TW' * Income / sum(Income);
    WelfarePctA2Ensemble2100h_TS = DQA2Ensemble2100h_TS' * Income / sum(Income);
    WelfarePctA1F12100h_TW = DQA1F12100h_TW' * Income / sum(Income);
    WelfarePctA1F12100h_TS = DQA1F12100h_TS' * Income / sum(Income);
    
    % Welfare impact of precip change
    if toggle_4<1       % precipitation not in model
        DQA2Ensemble2100h_P = zeros(N,1);
        DQA1F12100h_P = zeros(N,1);
    elseif toggle_4>=1 && toggle_4<=5    % avg annual precip
        DQA2Ensemble2100h_P = DeltaA2Ensemble2100_OW(:,1) * BetaHomog(S+KCG+1);
        DQA1F12100h_P = DeltaA1F12100_OW(:,1) * BetaHomog(S+KCG+1);
    else                            % summer and winter precip
        DQA2Ensemble2100h_P1 = DeltaA2Ensemble2100_OW3(:,1) * BetaHomog(S+KCG+1);
        DQA2Ensemble2100h_P2 = DeltaA2Ensemble2100_OW3(:,2) * BetaHomog(S+KCG+2);
        DQA2Ensemble2100h_P = DQA2Ensemble2100h_P1 + DQA2Ensemble2100h_P2;
        DQA1F12100h_P1 = DeltaA1F12100_OW3(:,1) * BetaHomog(S+KCG+1);
        DQA1F12100h_P2 = DeltaA1F12100_OW3(:,2) * BetaHomog(S+KCG+2);
        DQA1F12100h_P = DQA1F12100h_P1 + DQA1F12100h_P2;    end
    AggWelfareA2Ensemble2100h_P = DQA2Ensemble2100h_P' * Income;          % in billions
    WelfarePctA2Ensemble2100h_P = AggWelfareA2Ensemble2100h_P / sum(Income);
    AggWelfareA2Ensemble2100h_P = WelfarePctA2Ensemble2100h_P * USPersIncome2000;
    AggWelfareA1F12100h_P = DQA1F12100h_P' * Income;          % in billions
    WelfarePctA1F12100h_P = AggWelfareA1F12100h_P / sum(Income);
    AggWelfareA1F12100h_P = WelfarePctA1F12100h_P * USPersIncome2000;
    
    % Welfare impact of humid change
    if toggle_4<2        % humid not in model
        DQA2Ensemble2100h_H = zeros(N,1);
        DQA1F12100h_H = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=3    % average humid
        DQA2Ensemble2100h_H = DeltaA2Ensemble2100_OW(:,2) * BetaHomog(S+KCG+2);
        DQA1F12100h_H = DeltaA1F12100_OW(:,2) * BetaHomog(S+KCG+2);
    elseif toggle_4>=4 && toggle_4<=5    % summer humid (June-Aug)
        DQA2Ensemble2100h_H = DeltaA2Ensemble2100_OW2(:,2) * BetaHomog(S+KCG+2);
        DQA1F12100h_H = DeltaA1F12100_OW2(:,2) * BetaHomog(S+KCG+2);        
    else                            % summer and winter humid
        DQA2Ensemble2100h_H1 = DeltaA2Ensemble2100_OW3(:,3) * BetaHomog(S+KCG+3);
        DQA2Ensemble2100h_H2 = DeltaA2Ensemble2100_OW3(:,4) * BetaHomog(S+KCG+4);
        DQA2Ensemble2100h_H = DQA2Ensemble2100h_H1 + DQA2Ensemble2100h_H2;
        DQA1F12100h_H1 = DeltaA1F12100_OW3(:,3) * BetaHomog(S+KCG+3);
        DQA1F12100h_H2 = DeltaA1F12100_OW3(:,4) * BetaHomog(S+KCG+4);
        DQA1F12100h_H = DQA1F12100h_H1 + DQA1F12100h_H2;        
    end
    AggWelfareA2Ensemble2100h_H = DQA2Ensemble2100h_H' * Income;          % in billions
    WelfarePctA2Ensemble2100h_H = AggWelfareA2Ensemble2100h_H / sum(Income);
    AggWelfareA2Ensemble2100h_H = WelfarePctA2Ensemble2100h_H * USPersIncome2000;
    AggWelfareA1F12100h_H = DQA1F12100h_H' * Income;          % in billions
    WelfarePctA1F12100h_H = AggWelfareA1F12100h_H / sum(Income);
    AggWelfareA1F12100h_H = WelfarePctA1F12100h_H * USPersIncome2000;
    
    % Welfare impact of sunshine change
    % Be careful here. Control var in sunshine fraction but model change is
    % cloud fraction. So need to take negative.
    if toggle_4<2           % sunshine not in model
        DQA2Ensemble2100h_S = zeros(N,1);
        DQA1F12100h_S = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=5              % average annual sunshine
        DQA2Ensemble2100h_S = -DeltaA2Ensemble2100_OW(:,3) * BetaHomog(S+KCG+3);
        DQA1F12100h_S = -DeltaA1F12100_OW(:,3) * BetaHomog(S+KCG+3);
    else                            % summer and winter sunshine
        DQA2Ensemble2100h_S1 = -DeltaA2Ensemble2100_OW3(:,5) * BetaHomog(S+KCG+5);
        DQA2Ensemble2100h_S2 = -DeltaA2Ensemble2100_OW3(:,6) * BetaHomog(S+KCG+6);
        DQA2Ensemble2100h_S = DQA2Ensemble2100h_S1 + DQA2Ensemble2100h_S2;
        DQA1F12100h_S1 = -DeltaA1F12100_OW3(:,5) * BetaHomog(S+KCG+5);
        DQA1F12100h_S2 = -DeltaA1F12100_OW3(:,6) * BetaHomog(S+KCG+6);
        DQA1F12100h_S = DQA1F12100h_S1 + DQA1F12100h_S2;
    end
    AggWelfareA2Ensemble2100h_S = DQA2Ensemble2100h_S' * Income;          % in billions
    WelfarePctA2Ensemble2100h_S = AggWelfareA2Ensemble2100h_S / sum(Income);
    AggWelfareA2Ensemble2100h_S = WelfarePctA2Ensemble2100h_S * USPersIncome2000;
    AggWelfareA1F12100h_S = DQA1F12100h_S' * Income;          % in billions
    WelfarePctA1F12100h_S = AggWelfareA1F12100h_S / sum(Income);
    AggWelfareA1F12100h_S = WelfarePctA1F12100h_S * USPersIncome2000;  
    
    % Total "other weather" losses
    AggWelfareA2Ensemble2100h_OW = (DQA2Ensemble2100h_P+DQA2Ensemble2100h_H+DQA2Ensemble2100h_S)' * Income;          % in billions
    WelfarePctA2Ensemble2100h_OW = AggWelfareA2Ensemble2100h_OW / sum(Income);
    AggWelfareA2Ensemble2100h_OW = WelfarePctA2Ensemble2100h_OW * USPersIncome2000;
    AggWelfareA1F12100h_OW = (DQA1F12100h_P+DQA1F12100h_H+DQA1F12100h_S)' * Income;          % in billions
    WelfarePctA1F12100h_OW = AggWelfareA1F12100h_OW / sum(Income);
    AggWelfareA1F12100h_OW = WelfarePctA1F12100h_OW * USPersIncome2000;
    
    % Put together overall welfare impact and set up output
    DQA2Ensemble2100h = DQA2Ensemble2100h_T + DQA2Ensemble2100h_P + DQA2Ensemble2100h_H + DQA2Ensemble2100h_S;
    AggWelfareA2Ensemble2100h = DQA2Ensemble2100h' * Income;          % in billions
    WelfarePctA2Ensemble2100h = AggWelfareA2Ensemble2100h / sum(Income);
    AggWelfareA2Ensemble2100h = WelfarePctA2Ensemble2100h * USPersIncome2000;   
    DQA1F12100h = DQA1F12100h_T + DQA1F12100h_P + DQA1F12100h_H + DQA1F12100h_S;
    AggWelfareA1F12100h = DQA1F12100h' * Income;          % in billions
    WelfarePctA1F12100h = AggWelfareA1F12100h / sum(Income);
    AggWelfareA1F12100h = WelfarePctA1F12100h * USPersIncome2000;  
    PctLosingA2h = sum(PopStored.*(DQA2Ensemble2100h<0))/sum(PopStored);        % pct of population losing
    PctLosingA1F1h = sum(PopStored.*(DQA1F12100h<0))/sum(PopStored);        % pct of population losing
    
    
    % FINAL PIECE: MIGRATION
    ES = 8;      % Housing supply elasticity
    WelfarePctA2Ensemble2100hPop = DQA2Ensemble2100h' * PopStored / sum(PopStored); % pop-weighted
    DN = ES * (DQA2Ensemble2100h - WelfarePctA2Ensemble2100hPop);  % change in log population
    IncomeNew = exp(log(Income) + DN);
    WelfarePctA2Ensemble2100h_Migration = DQA2Ensemble2100h' * IncomeNew / sum(IncomeNew);
    PctPopChg = exp(DN) - 1;    % percent change in population

    % Store welfare estimates
    HOMOGOUTPUT(rep,:) = [WelfarePctA2Ensemble2100h WelfarePctA2Ensemble2100h_TW WelfarePctA2Ensemble2100h_TS WelfarePctA2Ensemble2100h_OW 0 PctLosingA2h WelfarePctA1F12100h       WelfarePctA1F12100h_TW       WelfarePctA1F12100h_TS       WelfarePctA1F12100h_OW       0 PctLosingA1F1h];
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % LOCAL LINEAR REGRESSION WITH BOOTSTRAP SAMPLE
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Set up vector of results
    BETA = zeros(N,KP);      % initialize

    
    % Take effects of linear parameters / coefficients out of the price
    % variable
    if FE==0 && DFE==0
        QRb = Qb - X(:,DropCols) * BetaHomog(DropCols);
    elseif FE==0 && DFE==1   
        QRb = Qb - [X(:,DropCols) DivFE] * [BetaHomog(DropCols); BetaHomog(K+1:KKK)];
    else
        QRb = Qb - [X(:,DropCols) StateFE] * [BetaHomog(DropCols); BetaHomog(K+1:KKK)];
    end

    % Population weight dep var 
    QRb = diag(Pop.^0.5) * QRb;
    
    % Now loop over each observation, doing a local regression at each
    % Weights don't change over bootstrap iterations since X's aren't being
    % re-drawn (though they still vary over the X's).
    for i = 1:N
        Z = zeros(N,KPW);     % initialize 
        % Need to build up kernel weighting matrix for all observations
        % if mod(i,500)==0; disp(i); end;       % show progress every 200 obs
        W = ones(N,1);     % initialize weights
        for k = 2:KPW            % ignore constant, so start at 2
            kk = KeepColsW(k);
            Z(:,k) = (XW(:,kk) - XW(i,kk)) / SigmaW(k);  % distance from current obs
        end
        Z = Z / h;
        Z = normpdf(Z,0,1);
        for k = 1:KPW
            W = W .* Z(:,k);
        end
        W = (W / h).^0.5;     % weights (must take square root given back-division below)

        % Now do weighted least squares    
        QHat = diag(W) * QRb;
        XHat = diag(W) * XR;    % XR already population weighted in original regression
        Qinv = inv(XHat' * XHat);
        
        BETA(i,:) = [Qinv * XHat' * QHat]';    % local coefs at location i
    end
    
    % RESTORE FULL SET OF BETA PARAMETERS (USING HOMOG PREF ESTIMATES)
    BETAALL = repmat(BetaHomog',N,1);       % N by KKK
    BETAALL(:,KeepCols) = BETA;             % Insert heterog prefs
    BETA = BETAALL;

    % Get value of an extra day at 65 degrees. Will need to subtract this off
    % of all estimates
    Val65 = BETA(:,1:S) * Basis(Bin65Ind,:)';

    % Get value of an extra day at each temperature in each location (rel to 65
    % degrees)
    BinValue = BETA(:,1:S) * Basis' - repmat(Val65,1,NbinInd);   % impact of one more day in each bin

        
    % If using splines, extrapolate MWTP at extremes
    if toggle_12>3 && toggle_12<50
        % Initialize valuations; BinValueh within present-day temp range, zero outside
        BinValueInd = BinValue;     % defined over present-day temp range
        BinValue = zeros(N,Nbin);   % defined over present & future temps (wider) 
        BinValue(:,Inds) = BinValueInd;
        if toggle_11==1      % do conf interval too
            TopIntInd = TopInt; BotIntInd = BotInt;
            TopInt = zeros(N,Nbin); BotInt = zeros(N,Nbin);
            TopInt(:,Inds) = TopIntInd; BotInt(:,Inds) = BotIntInd;
        end
        % Loop over Pumas, flattening valuations outside of the active
        % spline range
        for i = 1:N
            CenterS = CenterSh; LowS = LowSh; HighS = HighSh; MinS = MinSh; MaxS = MaxSh;
            ColdVal = BinValue(i,MinS); HotVal = BinValue(i,MaxS);    % valuations at bottom and top of temp dist
            BinValue(i,LowS) = ColdVal; BinValue(i,HighS) = HotVal;
            if toggle_11==1      % do conf interval too
                TopInt(i,LowS) = TopInt(i,MinS); TopInt(i,HighS) = TopInt(i,MaxS);
                BotInt(i,LowS) = BotInt(i,MinS); BotInt(i,HighS) = BotInt(i,MaxS);
            end    
        end
    end
        
    
    
    % Get weighted average MWTP across all PUMAs
    BinValueAvg = Income' * BinValue;
    BinValueAvg = BinValueAvg / sum(Income);
    BinValueAvg = BinValueAvg';
    
    
    
    % CONCAVITY CALCULATIONS
    SigmaX = std(XAll);    % std dev of each X variable (unweighted, including FE). For scaling   
    % Save the results of concavity of utility check in variable A
    Aij=zeros(N,N);     % Concavity needed to stop person at i from preferring j
    % Identify variables over which prefs are concave: all those that allow
    % random coefs (don't include constant term in linear spline models)
    if toggle_12<=3 || toggle_12>=50
        KeepColsC = KeepCols(2:KP);
    else
        KeepColsC = KeepCols(1:KP);     % cubic spline
    end
    % Get vector of unobserved attribute at each location
    Xi = zeros(N,1);
    for i = 1:N
        Xi(i) = Qb(i) - BETA(i,:) * XAll(i,:)';
    end
    % Loop over i to get Aij. Numerator is linear utility difference
    % between j and i. Denominator is quadratic characteristic difference
    for i = 1:N
        numerator = Qb(i) - Qb + Xi - Xi(i);
        numerator = numerator' + BETA(i,:) * XAll' - BETA(i,:) * XAll(i,:)';    % 1 X N vector
        sum_quadratic_term = sum(((ones(N,1)*XAll(i,KeepColsC)-XAll(:,KeepColsC)) ./ (ones(N,1)*SigmaX(KeepColsC))).^2, 2);
        sum_quadratic_term = sum_quadratic_term.^0.5;
        Aij(i,:) = numerator ./sum_quadratic_term';  % concavity needed to stop i from preferring j
        Aij(i,i) = 0; % Since denominator is 0 on the diagonals.
    end
    MaxConcav = max([Aij zeros(N,1)],[],2); % N by 1. Concavity needed at each i to stop i from preferring all other locations    
    
    
    
    
    
    
    
    
    %%%%%%%%%%%%
    % HETEROGENEOUS PREF WELFARE CALCULATIONS
    % Start with welfare impact of temperature change; linear terms
    DQA2Ensemble2100_T = zeros(N,1);     % initialize
    DQA1F12100_T = zeros(N,1);     % initialize
    for i = 1:N
        DQA2Ensemble2100_T(i) = BinValue(i,:) * DeltaA2Ensemble2100(i,:)';     % delta welfare as pct of income
        DQA1F12100_T(i) = BinValue(i,:) * DeltaA1F12100(i,:)';     % delta welfare as pct of income
    end
    
    % Additional Concavity Term. Do for each run of A2 scenario to account for
    % Jensen's inequality. All values are a percent of income
    if toggle_12>3 && toggle_12<50      % Flatten out extremes of spline and create new X's
        BasisF = zeros(Nbin,S);       % initalize full basis
        BasisF(Inds,:) = Basis;       % basis functions over present-day temperatures
        BasisF(LowS,:) = repmat(BasisF(MinS,:),length(LowS),1); % flatten basis at cold temps
        BasisF(HighS,:) = repmat(BasisF(MaxS,:),length(HighS),1); % flatten basis at cold temps
        XAllF(:,1:S) = Bin * BasisF;  % covariates corrected for flattened spline
        SigmaXF = std(XAllF);    % std dev of each X variable. For scaling   
        NonZeroXF = find(SigmaXF>0);        % covariates that are non-zero after flattening
        XAllF = XAllF(:,NonZeroXF); SigmaXF = SigmaXF(NonZeroXF);
        NonZeroSF = find(NonZeroXF<=S);       % elements of NonZeroX that pertain to splines
        BasisF = BasisF(:,NonZeroXF(NonZeroSF));  % non-degenerate basis functions
        SF = size(BasisF,2);            % number of remaining basis functions
    else
        SF = S; BasisF = Basis; XAllF = XAll; SigmaXF = SigmaX;
    end
    % Quadratic welfare losses
    % Include constant only if cubic spline model
    if toggle_12<=3 || toggle_12>=50; s=2; else; s=1; end;
    DQQA2Ensemble2100 = - sum(((BinA2Ensemble2100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
    DQQA1F12100 = - sum(((BinA1F12100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
%    DQQA2_a_2100 = - sum(((BinA2_a_2100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
%    DQQA2_b_2100 = - sum(((BinA2_b_2100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
%    DQQA2_c_2100 = - sum(((BinA2_c_2100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
%    DQQA2_d_2100 = - sum(((BinA2_d_2100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
%    DQQA2_e_2100 = - sum(((BinA2_e_2100*BasisF(:,s:SF) - XAllF(:,s:SF)) ./ (ones(N,1)*SigmaXF(s:SF))).^2,2).^0.5 .* MaxConcav;
    
    % Get change in marginal WTP given concavity and A2 ensemble warming
    BinValueQ = zeros(N,Nbin);      % initialize
    diff65 = (BinA2Ensemble2100-Bin)*BasisF.*repmat(BasisF(Bin65,:),N,1);   % N by SF. Change at 65F
    for b = 1:Nbin
        diff = (BinA2Ensemble2100-Bin)*BasisF.*repmat(BasisF(b,:),N,1);     % N by SF. Change at bin b
        BinValueQ(:,b) = -2*MaxConcav .* sum(((diff-diff65)./repmat(SigmaXF(1:SF).^2,N,1)),2);
    end
    clear diff diff65
        
    
    % Aggregate welfare in billions. Linear and quadratic terms
    AggWelfareA2Ensemble2100_T = DQA2Ensemble2100_T' * Income;          % in billions
    AggWelfareA1F12100_T = DQA1F12100_T' * Income;          % in billions
    AggWelfareQA2Ensemble2100 = DQQA2Ensemble2100' * Income;          % in billions
    AggWelfareQA1F12100 = DQQA1F12100' * Income;          % in billions
  %  AggWelfareQA2_a_2100 = DQQA2_a_2100' * Income;          % in billions
  %  AggWelfareQA2_b_2100 = DQQA2_b_2100' * Income;          % in billions
  %  AggWelfareQA2_c_2100 = DQQA2_c_2100' * Income;          % in billions
  %  AggWelfareQA2_d_2100 = DQQA2_d_2100' * Income;          % in billions
  %  AggWelfareQA2_e_2100 = DQQA2_e_2100' * Income;          % in billions
    
    % as a percent of total income. Linear and quadratic terms
    WelfarePctA2Ensemble2100_T = AggWelfareA2Ensemble2100_T / sum(Income);
    WelfarePctA1F12100_T = AggWelfareA1F12100_T / sum(Income);
    WelfarePctQA2Ensemble2100 = AggWelfareQA2Ensemble2100 / sum(Income);
    WelfarePctQA1F12100 = AggWelfareQA1F12100 / sum(Income);
 %   WelfarePctQA2_a_2100 = AggWelfareQA2_a_2100 / sum(Income);
 %   WelfarePctQA2_b_2100 = AggWelfareQA2_b_2100 / sum(Income);
 %   WelfarePctQA2_c_2100 = AggWelfareQA2_c_2100 / sum(Income);
 %   WelfarePctQA2_d_2100 = AggWelfareQA2_d_2100 / sum(Income);
 %   WelfarePctQA2_e_2100 = AggWelfareQA2_e_2100 / sum(Income);
    
    % Get welfare impact in $ using total US personal income
    AggWelfareA2Ensemble2100_T = WelfarePctA2Ensemble2100_T * USPersIncome2000;
    AggWelfareA1F12100_T = WelfarePctA1F12100_T * USPersIncome2000;
    AggWelfareQA2Ensemble2100 = WelfarePctQA2Ensemble2100 * USPersIncome2000;
    AggWelfareQA1F12100 = WelfarePctQA1F12100 * USPersIncome2000;
 %   AggWelfareQA2_a_2100 = WelfarePctQA2_a_2100 * USPersIncome2000;
 %   AggWelfareQA2_b_2100 = WelfarePctQA2_b_2100 * USPersIncome2000;
 %   AggWelfareQA2_c_2100 = WelfarePctQA2_c_2100 * USPersIncome2000;
 %   AggWelfareQA2_d_2100 = WelfarePctQA2_d_2100 * USPersIncome2000;
 %   AggWelfareQA2_e_2100 = WelfarePctQA2_e_2100 * USPersIncome2000;
    
    
    
    
    
    % Welfare impact of precip change
    if toggle_4<1       % precipitation not in model
        DQA2Ensemble2100_P = zeros(N,1);
        DQA1F12100_P = zeros(N,1);
    elseif toggle_4>=1 && toggle_4<=5    % avg annual precip
        DQA2Ensemble2100_P = DeltaA2Ensemble2100_OW(:,1) .* BETA(:,S+KCG+1);
        DQA1F12100_P = DeltaA1F12100_OW(:,1) .* BETA(:,S+KCG+1);
    else                            % summer and winter precip
        DQA2Ensemble2100_P1 = DeltaA2Ensemble2100_OW3(:,1) .* BETA(:,S+KCG+1);
        DQA2Ensemble2100_P2 = DeltaA2Ensemble2100_OW3(:,2) .* BETA(:,S+KCG+2);
        DQA2Ensemble2100_P = DQA2Ensemble2100h_P1 + DQA2Ensemble2100h_P2;
        DQA1F12100_P1 = DeltaA1F12100_OW3(:,1) .* BETA(:,S+KCG+1);
        DQA1F12100_P2 = DeltaA1F12100_OW3(:,2) .* BETA(:,S+KCG+2);
        DQA1F12100_P = DQA1F12100h_P1 + DQA1F12100h_P2;
    end
    AggWelfareA2Ensemble2100_P = DQA2Ensemble2100_P' * Income;          % in billions
    WelfarePctA2Ensemble2100_P = AggWelfareA2Ensemble2100_P / sum(Income);
    AggWelfareA2Ensemble2100_P = WelfarePctA2Ensemble2100_P * USPersIncome2000;
    AggWelfareA1F12100_P = DQA1F12100_P' * Income;          % in billions
    WelfarePctA1F12100_P = AggWelfareA1F12100_P / sum(Income);
    AggWelfareA1F12100_P = WelfarePctA1F12100_P * USPersIncome2000;
    
    % Welfare impact of summer humid change
    if toggle_4<2        % humid not in model
        DQA2Ensemble2100_H = zeros(N,1);
        DQA1F12100_H = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=3    % average annual humid
        DQA2Ensemble2100_H = DeltaA2Ensemble2100_OW(:,2) .* BETA(:,S+KCG+2);
        DQA1F12100_H = DeltaA1F12100_OW(:,2) .* BETA(:,S+KCG+2);
    elseif toggle_4>=4 && toggle_4<=5    % summer humid
        DQA2Ensemble2100_H = DeltaA2Ensemble2100_OW2(:,2) .* BETA(:,S+KCG+2);
        DQA1F12100_H = DeltaA1F12100_OW2(:,2) .* BETA(:,S+KCG+2);
    else                            % summer and winter humid
        DQA2Ensemble2100_H1 = DeltaA2Ensemble2100_OW3(:,3) .* BETA(:,S+KCG+3);
        DQA2Ensemble2100_H2 = DeltaA2Ensemble2100_OW3(:,4) .* BETA(:,S+KCG+4);
        DQA2Ensemble2100_H = DQA2Ensemble2100h_H1 + DQA2Ensemble2100h_H2;
        DQA1F12100_H1 = DeltaA1F12100_OW3(:,3) .* BETA(:,S+KCG+3);
        DQA1F12100_H2 = DeltaA1F12100_OW3(:,4) .* BETA(:,S+KCG+4);
        DQA1F12100_H = DQA1F12100h_H1 + DQA1F12100h_H2;
    end
    AggWelfareA2Ensemble2100_H = DQA2Ensemble2100_H' * Income;          % in billions
    WelfarePctA2Ensemble2100_H = AggWelfareA2Ensemble2100_H / sum(Income);
    AggWelfareA2Ensemble2100_H = WelfarePctA2Ensemble2100_H * USPersIncome2000;
    AggWelfareA1F12100_H = DQA1F12100_H' * Income;          % in billions
    WelfarePctA1F12100_H = AggWelfareA1F12100_H / sum(Income);
    AggWelfareA1F12100_H = WelfarePctA1F12100_H * USPersIncome2000;
    
    % Welfare impact of sunshine change
    % Be careful here. Control var in sunshine fraction but model change is
    % cloud fraction. So need to take negative.
    if toggle_4<2           % sunshine not in model
        DQA2Ensemble2100_S = zeros(N,1);
        DQA1F12100_S = zeros(N,1);
    elseif toggle_4>=2 && toggle_4<=5              % average annual sunshine
        DQA2Ensemble2100_S = -DeltaA2Ensemble2100_OW(:,3) .* BETA(:,S+KCG+3);
        DQA1F12100_S = -DeltaA1F12100_OW(:,3) .* BETA(:,S+KCG+3);
    else                            % summer and winter sunshine
        DQA2Ensemble2100_S1 = -DeltaA2Ensemble2100_OW3(:,5) .* BETA(:,S+KCG+5);
        DQA2Ensemble2100_S2 = -DeltaA2Ensemble2100_OW3(:,6) .* BETA(:,S+KCG+6);
        DQA2Ensemble2100_S = DQA2Ensemble2100h_S1 + DQA2Ensemble2100h_S2;
        DQA1F12100_S1 = -DeltaA1F12100_OW3(:,5) .* BETA(:,S+KCG+5);
        DQA1F12100_S2 = -DeltaA1F12100_OW3(:,6) .* BETA(:,S+KCG+6);
        DQA1F12100_S = DQA1F12100h_S1 + DQA1F12100h_S2;
    end
    AggWelfareA2Ensemble2100_S = DQA2Ensemble2100_S' * Income;          % in billions
    WelfarePctA2Ensemble2100_S = AggWelfareA2Ensemble2100_S / sum(Income);
    AggWelfareA2Ensemble2100_S = WelfarePctA2Ensemble2100_S * USPersIncome2000;
    AggWelfareA1F12100_S = DQA1F12100_S' * Income;          % in billions
    WelfarePctA1F12100_S = AggWelfareA1F12100_S / sum(Income);
    AggWelfareA1F12100_S = WelfarePctA1F12100_S * USPersIncome2000;
    
    % Total "other weather" losses
    AggWelfareA2Ensemble2100_OW = (DQA2Ensemble2100_P+DQA2Ensemble2100_H+DQA2Ensemble2100_S)' * Income;          % in billions
    WelfarePctA2Ensemble2100_OW = AggWelfareA2Ensemble2100_OW / sum(Income);
    AggWelfareA2Ensemble2100_OW = WelfarePctA2Ensemble2100_OW * USPersIncome2000;
    AggWelfareA1F12100_OW = (DQA1F12100_P+DQA1F12100_H+DQA1F12100_S)' * Income;          % in billions
    WelfarePctA1F12100_OW = AggWelfareA1F12100_OW / sum(Income);
    AggWelfareA1F12100_OW = WelfarePctA1F12100_OW * USPersIncome2000;
    
    % Put together overall welfare impact and set up output
    DQA2Ensemble2100 = DQA2Ensemble2100_T + DQA2Ensemble2100_P + DQA2Ensemble2100_H + DQA2Ensemble2100_S;
    AggWelfareA2Ensemble2100 = DQA2Ensemble2100' * Income;          % in billions
    WelfarePctA2Ensemble2100 = AggWelfareA2Ensemble2100 / sum(Income);
    AggWelfareA2Ensemble2100 = WelfarePctA2Ensemble2100 * USPersIncome2000;   
    DQA1F12100 = DQA1F12100_T + DQA1F12100_P + DQA1F12100_H + DQA1F12100_S;
    AggWelfareA1F12100 = DQA1F12100' * Income;          % in billions
    WelfarePctA1F12100 = AggWelfareA1F12100 / sum(Income);
    AggWelfareA1F12100 = WelfarePctA1F12100 * USPersIncome2000;  
    PctLosingA2 = sum(PopStored.*(DQA2Ensemble2100<0))/sum(PopStored);        % pct of population losing
    PctLosingA1F1 = sum(PopStored.*(DQA1F12100<0))/sum(PopStored);        % pct of population losing    
    
 
   
    % Store avg MWTPs and welfare estimates
    BinValueAvg_Bin40_365 = BinValueAvg(Bin40)*365; 
    BinValueAvg_Bin80_365 =  BinValueAvg(Bin80)*365;

    LINEAROUTPUT(rep,:) = [BinValueAvg_Bin40_365 BinValueAvg_Bin80_365 WelfarePctA2Ensemble2100 WelfarePctA1F12100];
    
    
    % Display repetition and output
    rep
    HOMOGOUTPUT(rep,:);
    LINEAROUTPUT(rep,:);
    
    end  
    
    
    
    
% BOOTSTRAP COMPLETE. WRAP UP AND STORE EVERYTHING IN TABLE-FRIENDLY FORMAT
    
std_HOMOGOUTPUT = std(HOMOGOUTPUT);
std_LINEAROUTPUT = std(LINEAROUTPUT);

OUTPUT1 = [HomogOutput(1) HomogOutput(3) HomogOutput(5:16); HomogOutput(2) HomogOutput(4)  std_HOMOGOUTPUT; LinearOutput 0 0 0 0 0 0 0 0 0 0 ; std_LINEAROUTPUT 0 0 0 0 0 0 0 0 0 0];

% Homogenous preference tables in paper
Tables4andA12 = OUTPUT1(1:2,:);
Tables4andA12 = reshape(Tables4andA12,28,1);
% Std. error of Total loss in 2008 billions can be directly computed by
% std. error of pct change: 
Tables4andA12(14) = Tables4andA12(6)*USPersIncome2000;
Tables4andA12(26) = Tables4andA12(18)*USPersIncome2000;

% Heterog prefs tables
Tables3and5 = OUTPUT1(3:4,1:4);
Tables3and5 = reshape(Tables3and5,8,1);

% csv output
dlmwrite('Tables4andA12_Q1_D2_FE0_S7_99.csv',Tables4andA12);
dlmwrite('Tables3and5_Q1_D2_FE0_S7_99.csv',Tables3and5);

save PartTWO_ALL_Q1_D2_FE0_S7_99.mat;
